Summary:
Analysis log for clinal fundulus project
grViz("digraph flowchart {
# node definitions with substituted label text
node [fontname = Helvetica, shape = rectangle]
tab1 [label = '@@1']
tab2 [label = '@@2']
tab3 [label = '@@3']
tab4 [label = '@@4']
tab5 [label = '@@5']
tab6 [label = '@@6']
tab7 [label = '@@7']
tab8 [label = '@@8']
# edge definitions with the node IDs
tab1 -> tab2 -> tab3 -> tab4 -> tab5 -> tab6 -> tab7
tab6 -> tab8
}
[1]: 'library prep (Elshire GBS)'
[2]: 'sequencing (illumina hiseq2000)'
[3]: 'qc checks and datawrangling (fastqc and bash)'
[4]: 'demultiplex and read filtering (stacks: process_radtags)'
[5]: 'read mapping (bwa_mem)'
[6]: 'genotype likelihoods (ANGSD)'
[7]: 'demography and population structure'
[8]: 'selection'
")
preserveb9b4482ed50a56cf
The clinal dataset is composed of 9 lanes of HiSeq2500 run in SE100. These reads include 1478 individuals from 5 separate experiments plus and additional set of lanes used for undersampled regions (CT, NC, SC and ME).
Lane names:
H7T6MADXX_1 H7T6MADXX_2 H9AREADXX_1 H9AREADXX_2 H9CWFADXX_1 H9WYGADXX_1 H9WYGADXX_2 HMMKLADXX_1 HMMKLADXX_2
Samples
sites<-read.table("./sites/sites.txt", header = T, sep = "\t")
kable(sites)
| Population_Full_Name | Abbreviation | N | Longitude | Latitude |
|---|---|---|---|---|
| Brayton Point, MA | BP | 50 | -71.18604 | 41.71250 |
| Clinton, CT | CT | 24 | -72.52448 | 41.27902 |
| Sapelo Isalnd Ferry Dock, Brunswick, GA | GA | 41 | -81.35452 | 31.44932 |
| Horseneck Beach, MA | HB | 50 | -71.02556 | 41.50449 |
| Kings Creek, Severn, VA | KC | 25 | -76.42462 | 37.29845 |
| Mattapoisett, MA | M | 31 | -70.81464 | 41.65875 |
| Mount Desert Island Bio Lab | MDIBL | 24 | -68.29441 | 44.42901 |
| Wiscassett, ME | ME | 27 | -69.66481 | 43.99695 |
| Mantoloking, NJ | MG | 331 | -74.07045 | 40.05046 |
| New Bedford Harbor | NBH, F, H, P, Syc | 150 | -70.90760 | 41.62486 |
| Manteo, NC | NC | 24 | -75.66737 | 35.90029 |
| Oyster Creek, NJ | OC | 47 | -74.18475 | 39.80866 |
| Rutgers, Tuckerton, NJ | RB | 423 | -74.32448 | 39.50878 |
| Charleston, SC | SC | 26 | -79.91624 | 32.75873 |
| Stone Harbor | SH | 62 | -74.76492 | 39.05666 |
| Slocum River, MA | SLC | 30 | -70.97109 | 41.52257 |
| Succotash Marsh, Matunuck, RI | SR | 49 | -71.52557 | 41.38235 |
| Wilmington, NC (Carolina Beach State Park) | WNC | 30 | -77.91932 | 34.04998 |
| grandis | grandis | 34 | -84.17740 | 30.07200 |
#!/bin/bash
# set the number of nodes
#SBATCH --nodes=1
# set the number of cpus
#SBATCH --cpus-per-task=6
# set max wall-clock time (D-HH:MM:SS)
#SBATCH --time=0-02:00:00
# set partition/queue to use
#SBATCH --partition=day-long-cpu
# set name of job
#SBATCH --job-name=fastqc
# set name of output file
#SBATCH --output=fastqc.out
# mail alert at start, end and abortion of execution
#SBATCH --mail-type=ALL
# send mail to this address
#SBATCH --mail-user=ddayan@clarku.edu
files="H7T6MADXX_1
H7T6MADXX_2
H9AREADXX_1
H9AREADXX_2
H9CWFADXX_1
H9WYGADXX_1
H9WYGADXX_2
HMMKLADXX_1
HMMKLADXX_2
"
for file in $files
do
/home/jgibbons/SOFTWARE/FastQC/fastqc -t 4 /home/ddayan/fundulus/seqdata/${file}_fastq.gz -o /home/ddayan/fundulus/seqdata/fastqc
done
fastqc looks good
We will use ANGSD to create population level allele frequency estimates rather than call indiviudal genotypes. The input files for ANGSD of indivudal level BAM files. Therefore the first step after QC check is to demultiplex, remove low quality reads and trim adapter and barcode sequences. Ths is accomplished with the process_radtags program from stacks v2.3.
Barcode keys (files used to demultiplex sequence data) are a bit complex given that some fish are sequenced twice and stacks requires each lane be run separately through process radtags. The solution is to give duplicate fish a unique ID suffix, but with a shared prefix, then cleaned demultiplexed reads from each fish in each lane are combined with their additional reads using prefix matching.
make new fish_ids, based on only population and library prep id, then mark duplicates (duplicate fish_ids occur because some individual fish are sequenced twice i.e. in different lanes)
barcodes<-read.table("./snp_calling/barcode_keys/cline_barcode_key.txt", sep = "\t", header = T)
tmp<-paste(barcodes$Pop , "_", barcodes$LibraryPrepID, sep = "")
barcodes$fish_id<-make.names(tmp, unique = TRUE)
write.csv(barcodes, row.names = FALSE, file = "./snp_calling/barcode_keys/cline_barcode_key.csv", quote = FALSE)
used python script on hand for slightly different format: create a unique single variable for sequencing lane (i.e. colbind flowcell and lane with a _) then run python script below
""" The input file has four columns, this script takes writes columns 2 and 3 (barcode and individual) to a new file based on the value of column 4."""
import csv
with open('./snp_calling/barcode_keys/cline_barcode_key.csv') as fin:
csvin = csv.DictReader(fin)
# Category -> open file lookup
outputs = {}
for row in csvin:
cat = row['Flowcell_Lane']
# Open a new file and write the header
if cat not in outputs:
fout = open('./snp_calling/barcode_keys/{}_key.csv'.format(cat), 'w')
dw = csv.DictWriter(fout, fieldnames=csvin.fieldnames)
dw.writeheader()
outputs[cat] = fout, dw
# Always write the row
outputs[cat][1].writerow(row)
# Close all the files
for fout, _ in outputs.values():
fout.close()
oops, only meant to keep barcode and sample id and need to write to tab delimited file
for i in ./snp_calling/barcode_keys/*csv
do
cut -d "," -f 2,10 $i > ${i%.csv}.tmp
done
for i in ./snp_calling/barcode_keys/*tmp
do
tr "," "\\t" < $i > ${i%.tmp}_barcodes.txt
done
clean up directory
rm ./snp_calling/barcode_keys/*.tmp
rm ./snp_calling/barcode_keys/*[12]_key.csv
remove first line from all files
sed -i -e "1d" ./meta_data/*_barcodes.txt
Summary
#!/bin/bash
# set max wall-clock time (D-HH:MM:SS)
#SBATCH --time=6-23:59:00
#SBATCH --cpus-per-task=38
# set partition/queue to use
#SBATCH --partition=week-long-cpu
# set name of job
#SBATCH --job-name=bwa_default
# set name of output file
#SBATCH --output=bwadefault.out
# send mail to this address
#SBATCH --mail-user=ddayan@clarku.edu
source /opt/samtools-1.6/source_me
for i in `find ../cleaned_tags/ -name '*.gz'`
do
/opt/bio-bwa/bwa mem -t 38 ../genome/f_het_genome_3.0.2.fa.gz ../cleaned/${i} | /opt/samtools-1.6/bin/samtools view -@ 16 -bSu - | /opt/samtools-1.6/bin/samtools sort -@ 16 - -o ./${i:0:-6}.$
done
What are the depths of mapped reads (from bam files)?
#!/bin/bash
# set max wall-clock time (D-HH:MM:SS)
#SBATCH --time=0-23:59:00
#SBATCH --cpus-per-task=1
# set partition/queue to use
#SBATCH --partition=day-long-cpu
# set name of job
#SBATCH --job-name=bam_depths
# set name of output file
#SBATCH --output=bam_depths
# send mail to this address
#SBATCH --mail-user=ddayan@clarku.edu
source /opt/samtools-1.6/source_me
samtools=/opt/samtools-1.6/bin/samtools
for i in *bam
do
$samtools depth $i | awk -v var="$i" '{sum += $3; n++} END {print " '${i%.bam}' ", sum / n}' >> depths.txt
done
bam_depths<-read.table("./snp_calling/QC/depths.txt", sep = " ")
bam_depths<-bam_depths[,-c(1,3)]
colnames(bam_depths)<-c("fish_id", "depth")
ggplot(data = bam_depths)+geom_density(aes(x = depth))+labs(title = "mean depth at mapped reads (bam) over all samples", x = "depth", y = "frequency")
#we could also easily split up the depth file here to look at population level variation in depth
let’s also look at overall number mapped reads per sample
#!/bin/bash
# set max wall-clock time (D-HH:MM:SS)
#SBATCH --time=0-23:59:00
#SBATCH --cpus-per-task=1
# set partition/queue to use
#SBATCH --partition=day-long-cpu
# set name of job
#SBATCH --job-name=bam_mapped_reads
# set name of output file
#SBATCH --output=bam_mapped_reads.out
source /opt/samtools-1.6/source_me
samtools=/opt/samtools-1.6/bin/samtools
for i in *bam
do
{ $samtools view -c -F 4 $i; echo ${i%.bam} ; } | tr "\n" " " >> mapped_reads.txt
done
mapped_reads<-read.table("./snp_calling/QC/mapped_reads.txt", sep = " ")
ggplot(data = mapped_reads)+geom_histogram(aes(x = log10(mapped_reads[,1])))+labs(title = " mapped reads in bam files", x = "log10 reads")
grandis
keeping grandis in the analysis may be useful to reconstruct the ancestral state of the genome and therefore create an unfolded SFS, but how well do the grandis reads align to the fundulus reference genome, let’s compare the bam files:
colnames(mapped_reads) <- c("reads", "sample")
mapped_reads$pop<-str_extract(mapped_reads$sample, "[^_]+")
mapped_reads %>%
group_by(pop) %>%
summarize(mean = mean(reads), sd = sd(reads), n = n())
ggplot()+geom_density(data = mapped_reads[(mapped_reads$pop!="grandis"),],aes(x=reads))+geom_density(data = mapped_reads[(mapped_reads$pop=="grandis"),],aes(x=reads), color = "red")+labs(title = "total number mapped reads", subtitle = "red - grandis ; black - heteroclitus")
this isn’t ideal though, lets try to figure out the PROPORTION of reads mapped to the reference genome
#!/bin/bash
# set max wall-clock time (D-HH:MM:SS)
#SBATCH --time=0-23:59:00
#SBATCH --cpus-per-task=1
# set partition/queue to use
#SBATCH --partition=day-long-cpu
# set name of job
#SBATCH --job-name=bam_mapped_reads
# set name of output file
#SBATCH --output=bam_mapped_reads.out
source /opt/samtools-1.6/source_me
samtools=/opt/samtools-1.6/bin/samtools
for i in *bam
do
{ $samtools view -c -F 0 $i; echo ${i%.bam} ; } | tr "\n" " " >> total_reads.txt
done
#total_reads <- read.table("./snp_calling/QC/total_reads.txt", sep = " ")
# make column of total reads, then proportion reads, plot again
ANGSD Analysis outline
We conduct two sets of ANGSD analyses. The first calculates genotype likelihoods. The second estimates the site frequency spectrum.
grViz("digraph flowchart {
# node definitions with substituted label text for placing code
node [fontname = Helvetica, shape = rectangle]
tab1 [label = '@@1']
tab2 [label = '@@2']
tab3 [label = '@@3']
tab4 [label = '@@4']
tab5 [label = '@@5']
tab6 [label = '@@6']
tab7 [label = '@@7']
tab8 [label = '@@8']
tab9 [label = '@@9']
# edge definitions with the node IDs
tab1 -> tab3;
tab1 -> tab2;
tab3 -> tab4 -> tab5
tab4 -> tab6
tab3 -> tab7
tab3 -> tab8
tab9 -> tab1
}
[1]: 'ANGSD -GL 1 -doMaf -doMajorMinor -glf 2'
[2]: 'Allele Frequencies'
[3]: 'beagle format GLs'
[4]: 'covariance matrix (PCangsd)'
[5]: 'PCA (R)'
[6]: 'RDA/other multivariate EAAs (R)'
[7]: 'LD (ngsLD)'
[8]: 'Admixture (ngsAdmix)'
[9]: 'bam files'
")
preservea5aff0014a782893
grViz("digraph flowchart {
# node definitions with substituted label text
node [fontname = Helvetica, shape = rectangle]
tab1 [label = '@@1']
tab2 [label = '@@2']
tab3 [label = '@@3']
tab4 [label = '@@4']
tab5 [label = '@@5']
tab6 [label = '@@6']
tab7 [label = '@@7']
tab8 [label = '@@8']
# edge definitions with the node IDs
tab1 -> tab2
tab2 -> tab3
tab3 -> tab4 -> tab5
tab3 -> tab6
tab3 -> tab7
tab8 -> tab1
}
[1]: 'ANGSD -GL 1 -doSAF unfolded'
[2]: 'SAF'
[3]: 'SFS (realSFS)'
[4]: 'Dadi conversion '
[5]: 'moments/Dadi'
[6]: 'theta, tajD'
[7]: 'FST'
[8]: 'bam files'
")
preservefcb942bf6cd370cc
Create lists of bams for all the separate populations
#in server directory containing bams
# for each population create a list of paths to bam files
find "$(cd ..; pwd)" -iname "BP*bam" > ../BP.bams
find "$(cd ..; pwd)" -iname "CT*bam" > ../CT.bams
find "$(cd ..; pwd)" -iname "GA*bam" > ../GA.bams
find "$(cd ..; pwd)" -iname "HB*bam" > ../HB.bams
find "$(cd ..; pwd)" -iname "KC*bam" > ../KC.bams
find "$(cd ..; pwd)" -iname "M_*bam" > ../M.bams
find "$(cd ..; pwd)" -iname "MDIBL_*bam" > ../MDIBL.bams
find "$(cd ..; pwd)" -iname "ME_*bam" > ../ME.bams
find "$(cd ..; pwd)" -iname "MG_*bam" > ../MG.bams
find "$(cd ..; pwd)" -iname "NBH_*bam" > ../NBH.bams
find "$(cd ..; pwd)" -iname "F_*bam" > ../F.bams
find "$(cd ..; pwd)" -iname "H_*bam" > ../H.bams
find "$(cd ..; pwd)" -iname "P_*bam" > ../P.bams
find "$(cd ..; pwd)" -iname "SYC_*bam" > ../SYC.bams
find "$(cd ..; pwd)" -iname "NC_*bam" > ../NC.bams
find "$(cd ..; pwd)" -iname "OC_*bam" > ../OC.bams
find "$(cd ..; pwd)" -iname "RB_*bam" > ../RB.bams
find "$(cd ..; pwd)" -iname "SC_*bam" > ../SC.bams
find "$(cd ..; pwd)" -iname "SH_*bam" > ../SH.bams
find "$(cd ..; pwd)" -iname "SLC_*bam" > ../SLC.bams
find "$(cd ..; pwd)" -iname "SR_*bam" > ../SR.bams
find "$(cd ..; pwd)" -iname "WNC_*bam" > ../WNC.bams
find "$(cd ..; pwd)" -iname "grandis_*bam" > ../grandis.bams
cat *bams > ALL.bamlist
first compress the genome and create a compression index
#!/bin/bash
# set max wall-clock time (D-HH:MM:SS)
#SBATCH --time=0-23:59:00
#SBATCH --cpus-per-task=10
# set partition/queue to use
#SBATCH --partition=day-long-cpu
# set name of job
#SBATCH --job-name=bgzip_index
# set name of output file
#SBATCH --output=bgzip.out
# send mail to this address
#SBATCH --mail-user=ddayan@clarku.edu
source /opt/samtools-1.6/source_me
/home/ddayan/software/htslib-1.9/bgzip -i -@ 10 f_het_genome_3.0.2.fa
then create a samtools index
samtools faidx f_het_genome_3.0.2.fa.gz
also need to index all the bam files
#!/bin/bash
# set max wall-clock time (D-HH:MM:SS)
#SBATCH --time=0-23:59:00
#SBATCH --cpus-per-task=10
# set partition/queue to use
#SBATCH --partition=day-long-cpu
# set name of job
#SBATCH --job-name=bgzip_index
# set name of output file
#SBATCH --output=bgzip.out
# send mail to this address
#SBATCH --mail-user=ddayan@clarku.edu
source /opt/samtools-1.6/source_me
samtools=/opt/samtools-1.6/bin/samtools
for i in ./*.bam;do $samtools index $i;done
First lets take a look at coverage and allele freqeuncy to get an idea what the data looks like prior to filtering
rationale for options:
command options
-P 10: use 10 threads
-ref: path to reference genome
-out: output
-r : run on only one contig (to check speed)
filtering
-uniqueOnly: use only uniquley mapping reads
-remove_bads: discard bad reads
-trim 0: perform no trimming
-C exclude reads with excessive mismatches
-baq avoid SNPs called due to indels
-minmapQ exclude reads with poor mapping
-maxDepth - r)ead depth above this level are binned (set high to llok at actual depth distrubution
dos
doqsdist calulates distribution of quality scores of used reads
dodepth calculates distrubtion of read depths
docounts needed for depth calcuation
#!/bin/bash
# set max wall-clock time (D-HH:MM:SS)
#SBATCH --time=6-23:59:00
#SBATCH --cpus-per-task=10
# set partition/queue to use
#SBATCH --partition=week-long-cpu
# set name of job
#SBATCH --job-name=bwa_default
# set name of output file
#SBATCH --output=bwadefault.out
# send mail to this address
#SBATCH --mail-user=ddayan@clarku.edu
ANGSD="/opt/angsd0.928/angsd/angsd"
REF="/home/ddayan/fundulus/genome/f_het_genome_3.0.2.fa.gz"
$ANGSD -P 10 -b ../meta_data/ALL.bamlist -ref $REF -out ./Results/ALL.qc \
-uniqueOnly 1 -remove_bads 1 -trim 0 -C 50 -baq 1 -maxDepth 1000 \
-minMapQ 20 -doQsDist 1 -doDepth 1 -doCounts 1
results of initial ANGSD run Q scores and global depth
## barplot q-scores
qs <- read.table(file = "./snp_calling/ANGSD_outputs/qc_check/ALL.qc.qs", head=T, stringsAsFactors=F)
barplot(height=qs$counts, names.arg=qs$qscore, xlab="Q-score", ylab="Counts", main = "Q-score distribution")
## global depth
dep <- as.numeric(scan(file = "./snp_calling/ANGSD_outputs/qc_check/ALL.qc.depthGlobal",what="char", quiet=T))
dep_freq<-as.data.frame(cbind(dep, seq(1:length(dep))))
colnames(dep_freq)<-c("count", "depth")
dep_freq$depth<-dep_freq$depth-1
ggplot(data = dep_freq, aes(depth, count))+geom_bar(stat = "identity")+labs(x = "Depth", title = "global depths (summed across all fish)")
ggplot(data = dep_freq, aes(depth, count))+geom_bar(stat = "identity")+labs(x = "Depth", title = "subset of plot above")+xlim(0, 100)
Depths for a handful of individuals:
## individual depth
idep<-read.table("./snp_calling/ANGSD_outputs/qc_check/ALL.qc.depthSample")
#look at read depth histogram for the first tenindividual
i1<-as.data.frame(t(idep[1,]))
i1$depth<-c(0:999, by = 1)
colnames(i1)<-c("log10count", "depth")
i1$log10count <- log10(i1$log10count)
ggplot(data = i1, aes(depth, log10count))+geom_bar(stat = "identity")+labs(x = "Depth")
i200<-as.data.frame(t(idep[1,]))
i200$depth<-c(0:999, by = 1)
colnames(i200)<-c("log10count", "depth")
i200$log10count <- log10(i200$log10count)
ggplot(data = i200, aes(depth, log10count))+geom_bar(stat = "identity")+labs(x = "Depth")
i505<-as.data.frame(t(idep[1,]))
i505$depth<-c(0:999, by = 1)
colnames(i505)<-c("log10count", "depth")
i505$log10count <- log10(i505$log10count)
ggplot(data = i200, aes(depth, log10count))+geom_bar(stat = "identity")+labs(x = "Depth", title = )
other salient results Total number of sites analyzed: 567547 Number of sites retained after filtering: 559561
thoughts here: - most sites have very limited depth, with 1480 samples, the global depth is largely less than 1000
- 59000 sites (about 10%) have more than 1000 reads, are these driven by a moderate depth in a lot of individuals or extremely high depth in just a few?
- lets run this again but only retain sites present in 80% of individuals
#!/bin/bash
# set max wall-clock time (D-HH:MM:SS)
#SBATCH --time=6-23:59:00
#SBATCH --cpus-per-task=10
# set partition/queue to use
#SBATCH --partition=week-long-cpu
# set name of job
#SBATCH --job-name=bwa_default
# set name of output file
#SBATCH --output=bwadefault.out
# send mail to this address
#SBATCH --mail-user=ddayan@clarku.edu
ANGSD="/opt/angsd0.928/angsd/angsd"
REF="/home/ddayan/fundulus/genome/f_het_genome_3.0.2.fa.gz"
$ANGSD -P 10 -b ../meta_data/ALL.bamlist -ref $REF -out ./Results/ALL.qc -r KN811289.1 \
-uniqueOnly 1 -remove_bads 1 -trim 0 -C 50 -baq 1 -maxDepth 1000 -minInd 1184 \
-minMapQ 20 -doQsDist 1 -doDepth 1 -doCounts 1
no sites were retained with this level of filtering, lets scale up to the ful genome and see how many sites with coverage in 80% of individuals we get: 2702 sites …
this ANGSD output is very sparse, how many loci are these sites spread across? how many reads are they based on before and after filtering? what are the mapping qualities? why is the depth data presented in the least convenient imaginable format? I miss stacks and tassel…
lets try again but a little less stringent, here we retain sites with coverage in 70% of fundulus heterociltis in the dataset (1010),
#!/bin/bash
# set max wall-clock time (D-HH:MM:SS)
#SBATCH --time=0-23:59:00
#SBATCH --cpus-per-task=10
# set partition/queue to use
#SBATCH --partition=day-long-cpu
# set name of job
#SBATCH --job-name=angsd_0.7
# set name of output file
#SBATCH --output=angsd_0.7.out
ANGSD="/opt/angsd0.928/angsd/angsd"
REF="/home/ddayan/fundulus/genome/f_het_genome_3.0.2.fa.gz"
$ANGSD -P 10 -b ../meta_data/ALL.bamlist -ref $REF -out ./Results/ALL.qc -r KN811289.1 \
-uniqueOnly 1 -remove_bads 1 -trim 0 -C 50 -baq 1 -minInd 1010 \
-minMapQ 20
how many loci is this over? we can use gatk countloci for this
prep:
GATK="/opt/java/jdk1.8.0_121/bin/java -jar /opt/Gatk/GenomeAnalysisTK.jar"
#create dictionary for gatk
java -jar /opt/picard-tools-1.119/CreateSequenceDictionary.jar R= fundulus/genome/f_het_genome_3.0.2.fa O= fundulus/genome/f_het_genome_3.0.2.dict
#create an index for the uncompressed genome file (gatk only run on decompressed reference genomes)
samtools faidx f_het_genome_3.0.2.fa
#now we can finally try to run GATK
$GATK -T CountLoci -R f_het_genome_3.0.2.fa -I ../aligned/OC_421.merged.bam
####agh this doesn't work because read groups were not retained in the bam files...
Here we’re going to finally run ANGSD
#!/bin/bash
# set max wall-clock time (D-HH:MM:SS)
#SBATCH --time=6-23:59:00
#SBATCH --cpus-per-task=38
# set partition/queue to use
#SBATCH --partition=week-long-cpu
# set name of job
#SBATCH --job-name=angsd_GL_0.1
# set name of output file
#SBATCH --output=angsd_GL_0.1.out
ANGSD="/opt/angsd0.928/angsd/angsd"
REF="/home/ddayan/fundulus/genome/f_het_genome_3.0.2.fa.gz"
FILTERS="-minInd 740 -uniqueOnly 1 -remove_bads 1 -trim 0 -C 50 -baq 1 -skipTriallelic 1 -SNP_pval 1e-6 -minMapQ 20 -maxDepth 14800"
DOS="-GL 1 -doMajorMinor 1 -doMaf 2 -doGlf 2 -doCounts 1 -doDepth 1"
$ANGSD -P 38 -b ../meta_data/ALL.bamlist -ref $REF -out ./Results/ALL.qc \
$FILTERS $DOS \
Total number of sites analyzed: 80148766
Number of sites retained after filtering: 62571
Coverage
## global depth
dep <- as.numeric(scan(file = "./snp_calling/ANGSD_outputs/run_0.1//ALL.qc.depthGlobal",what="char", quiet=T))
dep_freq<-as.data.frame(cbind(dep, seq(1:length(dep))))
colnames(dep_freq)<-c("count", "depth")
dep_freq$depth<-dep_freq$depth-1
dep_freq$log_count<-log10(dep_freq$count)
ggplot(data = dep_freq, aes(depth, log_count))+geom_bar(stat = "identity")+labs(x = "Depth", title = "global depths (summed across all fish)")
Of the 62571 sites, 47984 (76%) have more than 14800 reads and show up as a single thin bar at 14800 in this plot
## individual depth
idep<-read.table("./snp_calling/ANGSD_outputs/run_0.1/ALL.qc.depthSample")
#look at read depth histogram for the first tenindividual
#i100<-as.data.frame(t(idep[1,]))
#i100$depth<-c(0:14799, by = 1)
#colnames(i100)<-c("log10count", "depth")
#i100$log10count <- log10(i100$log10count)
#ggplot(data = i100, aes(depth, log10count))+geom_bar(stat = "identity")+labs(x = "Depth")+xlim(c(0,500))
i100<-as.data.frame(t(idep[100,]))
i100$depth<-c(0:14799, by = 1)
i100$log10count <- log10(i100$`1`)
colnames(i100)<-c("count", "depth", "log10count")
#ggplot(data = i100, aes(depth, log10count))+geom_bar(stat = "identity")+labs(x = "Depth")+xlim(c(0,500))
ggplot(data = i100, aes(depth, count))+geom_bar(stat = "identity")+labs(x = "Depth")+xlim(c(0,400))+ylim(c(0,4000))+labs(title = "individual 100")
## Warning: Removed 14400 rows containing missing values (position_stack).
## Warning: Removed 1 rows containing missing values (geom_bar).
i1000<-as.data.frame(t(idep[1000,]))
i1000$depth<-c(0:14799, by = 1)
i1000$log10count <- log10(i1000$`1000`)
colnames(i1000)<-c("count", "depth", "log10count")
#ggplot(data = i1000, aes(depth, log10count))+geom_bar(stat = "identity")+labs(x = "Depth")+xlim(c(0,1000))
ggplot(data = i1000, aes(depth, count))+geom_bar(stat = "identity")+labs(x = "Depth")+xlim(c(0,400))+ylim(c(0,2000))+labs(title = "individual 1000")
## Warning: Removed 14400 rows containing missing values (position_stack).
## Warning: Removed 1 rows containing missing values (geom_bar).
i500<-as.data.frame(t(idep[500,]))
i500$depth<-c(0:14799, by = 1)
i500$log10count <- log10(i500$`500`)
colnames(i500)<-c("count", "depth", "log10count")
#ggplot(data = i500, aes(depth, log10count))+geom_bar(stat = "identity")+labs(x = "Depth")+xlim(c(0,500))
ggplot(data = i500, aes(depth, count))+geom_bar(stat = "identity")+labs(x = "Depth")+xlim(c(0,300))+ylim(c(0,3000))+labs(title = "individual 500")
## Warning: Removed 14500 rows containing missing values (position_stack).
## Warning: Removed 1 rows containing missing values (geom_bar).
i148<-as.data.frame(t(idep[148,]))
i148$depth<-c(0:14799, by = 1)
i148$log10count <- log10(i148$`148`)
colnames(i148)<-c("count", "depth", "log10count")
#ggplot(data = i148, aes(depth, log10count))+geom_bar(stat = "identity")+labs(x = "Depth")+xlim(c(0,148))
ggplot(data = i148, aes(depth, count))+geom_bar(stat = "identity")+labs(x = "Depth")+xlim(c(0,300))+ylim(c(0,3000))+labs(title = "grandis individual")
## Warning: Removed 14500 rows containing missing values (position_stack).
## Warning: Removed 1 rows containing missing values (geom_bar).
if we look at a handful of depth for individual fish (reads per locus per individual) we observe a nice distribution with mostly >20x coverage
a note on installation (code chunk below):
#didn't want to wait for HPC to install pcangsd so did it locally. the problem (as always) is that pcangsd has depnedinecies that get install to local python but slurm calls a different python. used the tool virtualenv to get around this
#create the virtual environment in the directory containing pcagnsd
virtualenv software/pcangsd/
cd software/pcangsd/
#ctivate it
source bin/activate
#install dependcies and angsd
pip install cython
pip install numpy
python setup.py build_ext --inplace
pip install -r requirements.txt
#check install
python pcangsd.py
#deactive virtualenv
#when running in batch script make sure to source this virtual env before executing the actual command
deactivate
ok now lets try a PCA
#!/bin/bash
# set max wall-clock time (D-HH:MM:SS)
#SBATCH --time=6-23:59:00
#SBATCH --cpus-per-task=38
# set partition/queue to use
#SBATCH --partition=week-long-cpu
# set name of job
#SBATCH --job-name=pcangsd_0.1
# set name of output file
#SBATCH --output=pcangsd_0.1.out
source /home/ddayan/software/pcangsd/bin/activate
python /home/ddayan/software/pcangsd/pcangsd.py -beagle ../Results/ALL_gl_0.1.beagle.gz -o pcangsd.gl_0.1 -threads 38
Let’s looks at these PCA results
cov_mat<-read.table("./pcangsd/pcangsd.gl_0.1.cov", sep = " ")
filenames = scan("./snp_calling/ANGSD_outputs/ALL.bamlist", what="character")
fish_names = gsub("/home/ddayan/fundulus/aligned/", "", gsub(".merged", "", gsub(".bam", "", filenames)))
pc1<-prcomp(cov_mat)
pc_inds<-as.data.frame(pc1$x)
row.names(pc_inds)<-fish_names
pc_inds$pop<-gsub("_\\d+", "", fish_names)
#caution this deletes al nbh pops except nbh, fix later
pc_inds <- merge(pc_inds, sites[,c("Abbreviation", "Latitude")], by.x = "pop" , by.y ="Abbreviation")
ggplot(data = pc_inds) + geom_point(aes(x = PC1, y = PC2, color = Latitude))+scale_color_gradientn(colours = rainbow(5))#+stat_ellipse(aes(x = PC1, y = PC2, group = pop),type = "norm")
plot_ly(x=pc_inds$PC1, y=pc_inds$PC2, z=pc_inds$PC3, type="scatter3d", mode="markers", color=pc_inds$pop)
preserved7105fe5b51186fa
hover mouse over points for population id - population codes are as in the table “sites” under the section “data description” at top of document
it looks like there are some QC issues to take on (or a grandis individual is more similar to f het than other grandis), but that this approach is at least working well enough to differentiate pops in PC space
grandis
we will use the grandis reads to create a grandis consensus sequence (dofasta in angsd) which will be used as the ancestral state for unfolded SFS estimation, but we will not estimate genotype likelihoods in grandis as this will influence SNP probabilities and it seems that including grandis dominates the structure found by PCA
HWE
no HWE filtering, filtering out paralogues by using excess coverage instead - rationale here is that HWE is likely to be skewed given the amount of population strcutre in the sample and that filtering by HWE stratified by population is not possible in the ANGSD pipeline
MAF
MAF filtering will skew the SFS and should not be applied to the methods that rely on SFS (demography inference, FST), but the inclusion of singletons called genotyped datasets can confound STRUCTURE (little effect on multivariate population structure algorithms e.g. PCA DAPC). My assumption is that the use of genotype likelihoods instead of called genotypes should reduce the impact of singletons or low frequency alleles (erroneous or true) on the inference of structure using model based approaches like STRUCTURE, and the inclusion of rare alleles will improve discrimination among recently diverged populations (or populations with rare gene flow) when using multivariate approaches
max_coverage
will set max coverage at 2.5x the modal coverage per locus, this means will have to run ANGSD again with a higher maxdepth cutoff for the coverage output in order to estimate it…
min_coverage_locus
no hard minimum coverage cutoff, instead use stringent snp_pval cutoff and minimum number of individuals (we require that a site be present in 50% of samples - 686)
minimum_coverage_ind
the mean coverage for the individuals (at the bam level) is ~730k, but a handful of individuals have very low coverege (see bam stats section), will remove these individuals BEFORE ANGSD analysis by editing bam_lists, removed individuals in the lowest 5% of the coverage distribution
snp_pval
impose stringent probablility cutoff for reatining a SNP: 1e-6
quality_and_mapping
-uniqueOnly 1 -remove_bads 1 -C 50 -baq 1 -skipTriallelic 1 -minMapQ 20
version history 1.0 first run 1.1 removed outlier high coverage 1.2 start from scratch with mislabeled GA fish removed - GA_3318 clustered with grandis in PC space and also was it’s own group in PCA run on genotype likelihoods v1.1 1.3 coverage filtered output of 1.2
Make list of bam files, filtering out grandis and individuals in the lowest 5% of reads:
#get rid of individuals with low sequencing
bad_inds <- mapped_reads[mapped_reads$reads<quantile(mapped_reads$reads, 0.05),2]
#also get rid of mislabeled individual - pca_1.1 showed that "GA_3318" is most likely a grandis
bad_inds<-append(bad_inds, mapped_reads[140,2])
bad_inds <- mapped_reads[bad_inds,2]
#write.table(bad_inds, file = "./snp_calling/QC/bad_inds.txt", quote = FALSE, row.names = FALSE)
grep -v -f ./snp_calling/QC/bad_inds.txt ./snp_calling/ANGSD_outputs/ALL.bamlist | grep -v 'grandis' > ./snp_calling/QC/good_no_grandis_2.bamlist
This removed all grandis, 72 individuals with too few reads and 1 individual labeled as GA but actually grandis, leaving a final dataset of 1371 individuals.
Make genotype likelihood estimation v1.2
#!/bin/bash
# set max wall-clock time (D-HH:MM:SS)
#SBATCH --time=6-23:59:00
#SBATCH --cpus-per-task=38
# set partition/queue to use
#SBATCH --partition=week-long-cpu
# set name of job
#SBATCH --job-name=angsd_GL_1.0
# set name of output file
#SBATCH --output=angsd_GL_1.0.out
ANGSD="/opt/angsd0.928/angsd/angsd"
REF="/home/ddayan/fundulus/genome/f_het_genome_3.0.2.fa.gz"
FILTERS="-minInd 686 -uniqueOnly 1 -remove_bads 1 -trim 0 -C 50 -baq 1 -skipTriallelic 1 -SNP_pval 1e-6 -minMapQ 20"
DOS="-GL 1 -doMajorMinor 1 -doMaf 2 -doGlf 2 -doCounts 1 -dumpCounts 2"
$ANGSD -P 38 -b ../meta_data/good_no_grandis_2.bamlist -ref $REF -out ./Results/gl_1.2 \
$FILTERS $DOS \
Run 1.0: Total number of sites analyzed: 78663923 Number of sites retained after filtering: 77717
Run 1.2 Total number of sites analyzed: 78419827 Number of sites retained after filtering: 75430
Extract and plot coverage data
#check the format before running this
depths_loci<-read.table("./snp_calling/QC/gl_1.2.pos", header = T, sep = "\t")
#mode calculation
getMode <- function(x) {
keys <- unique(x)
keys[which.max(tabulate(match(x, keys)))]
}
modal_depth <- getMode(depths_loci$totDepth)
ggplot(data = depths_loci) + geom_histogram(aes(x = totDepth))+xlim(c(0,75000)) + geom_vline(aes(xintercept = modal_depth), color = "red") + geom_vline(aes(xintercept = 2.5*modal_depth), color = "blue")
The modal depth per locus is modal_depth. All loci with 2.5x greater read depth are then filtered out.
#note, remove comments if need to write the outputs to file
#find loci with >2.5x the mode and filter then write out to loci list which will be used to filter the GL outputs (filtered datasets are named 1.1)
sites_to_keep<-subset(depths_loci, totDepth < 2.5*modal_depth)
#write.table(sites_to_keep[,c(1,2)], file = "./snp_calling/QC/sites_to_keep.txt", row.names = FALSE, quote = FALSE, col.names = FALSE)
#get bad sites and store in beagle call format
sites_to_remove<-subset(depths_loci, totDepth > 2.5*modal_depth)
sites_to_remove$beagle_id<-paste(sites_to_remove$chr, "_", sites_to_remove$pos, sep = "")
#write.table(sites_to_remove[,c(4)], file = "./snp_calling/QC/sites_to_remove.txt", row.names = FALSE, quote = FALSE, col.names = FALSE)
GL_1.2 had nrow(depths_loci) SNPs, after fitlering out SNPs with greater than 2.5x the mode coverage (modal_depthx = mode) there are nrow(sites_to_keep) SNPs
#filtering the datasets
zcat gl_1.2.beagle.gz > gl_1.2.beagle
awk 'FNR==NR { a[$1]; next } !($1 in a){print $0}' sites_to_remove.txt gl_1.2.beagle > gl_1.3.beagle
bgzip gl_1.3.beagle
#!/bin/bash
# set max wall-clock time (D-HH:MM:SS)
#SBATCH --time=6-23:59:00
#SBATCH --cpus-per-task=38
# set partition/queue to use
#SBATCH --partition=week-long-cpu
# set name of job
#SBATCH --job-name=pcangsd_1.3
# set name of output file
#SBATCH --output=pcangsd_1.3.out
source /home/ddayan/software/pcangsd/bin/activate
python /home/ddayan/software/pcangsd/pcangsd.py -beagle ../Results/gl_1.3.beagle.gz -o pcangsd.gl_1.1 -threads 38 -minMaf 0.0
PCANGSD run info:
Read 1371 samples and 73375 sites
Estimating population allele frequencies EM (MAF) converged at iteration: 12
Estimating covariance matrix Using 5 principal components (MAP test)
cov_mat<-read.table("./pcangsd/pcangsd.gl_1.3.cov", sep = " ")
filenames <- mapped_reads[! mapped_reads$sample %in% bad_inds,]
filenames <- filenames[!grepl("grandis", filenames$sample),]
filenames <- filenames[-c(1372, 1373),]
fish_names = gsub(".merged", "", filenames$sample)
pc1<-prcomp(cov_mat)
save(pc1, file = "./pcangsd/pca.R")
pc_inds<-as.data.frame(pc1$x)
row.names(pc_inds)<-fish_names
pc_inds$pop<-gsub("_\\d+", "", fish_names)
# write.table(pc_inds, file = "./pcangsd/pca.txt", quote = FALSE) #commented to avoid overwriting, run if doing first run through analysis
#caution this deletes al nbh pops except nbh, fix later
pc_inds <- merge(pc_inds, sites[,c("Abbreviation", "Latitude")], by.x = "pop" , by.y ="Abbreviation")
ggplot(data = pc_inds) + geom_point(aes(x = PC1, y = PC2, color = Latitude))+scale_color_gradientn(colours = rainbow(5))#+stat_ellipse(aes(x = PC1, y = PC2, group = pop),type = "norm")
plot_ly(x=pc_inds$PC1, y=pc_inds$PC2, z=pc_inds$PC3, type="scatter3d", mode="markers", color=pc_inds$pop)
preserve506392caba688e5e
There is 1 outlier fish (“GA_3318”). Given that this fish was caught in Georgia, it is likely a misidentified grandis is needs to be excluded from the analysis. New GL file and all resulting analysis labeled 1.2
NGSadmix is used to estimate ancestry for K = 1-6 with 10 replicate runs and best k determined by evanno
#!/bin/bash
# set max wall-clock time (D-HH:MM:SS)
#SBATCH --time=6-23:59:00
#SBATCH --cpus-per-task=10
# set partition/queue to use
#SBATCH --partition=week-long-cpu
# set name of job
#SBATCH --job-name=ngsadmix_1.3_k9
# set name of output file
#SBATCH --output=ngsadmix_1.3_k9
NGSadmix="/opt/angsd0.928/angsd/misc/NGSadmix"
BEAGLE="/home/ddayan/fundulus/ANGSD/Results/gl_1.3.beagle.gz"
Out_Prefix="1.3_k9"
K="9"
$NGSadmix -likes $BEAGLE -K $K -o $Out_Prefix.1 -P 10 -minMaf 0
$NGSadmix -likes $BEAGLE -K $K -o $Out_Prefix.2 -P 10 -minMaf 0
$NGSadmix -likes $BEAGLE -K $K -o $Out_Prefix.3 -P 10 -minMaf 0
$NGSadmix -likes $BEAGLE -K $K -o $Out_Prefix.4 -P 10 -minMaf 0
$NGSadmix -likes $BEAGLE -K $K -o $Out_Prefix.5 -P 10 -minMaf 0
$NGSadmix -likes $BEAGLE -K $K -o $Out_Prefix.6 -P 10 -minMaf 0
$NGSadmix -likes $BEAGLE -K $K -o $Out_Prefix.7 -P 10 -minMaf 0
$NGSadmix -likes $BEAGLE -K $K -o $Out_Prefix.8 -P 10 -minMaf 0
$NGSadmix -likes $BEAGLE -K $K -o $Out_Prefix.9 -P 10 -minMaf 0
$NGSadmix -likes $BEAGLE -K $K -o $Out_Prefix.10 -P 10 -minMaf 0
Plot admixture portions for a single run of k = 2
#make popfile
pops<- as.data.frame(cbind(fish_names, gsub("_\\d+", "", fish_names)))
colnames(pops) <- c("sample", "pop")
#plot a admixture plot
q<-read.table("./ngs_admix/sandbox/outFileName.qopt")
#set order of poulations by latitude
q$sample <- pops$sample
q$pop <- pops$pop
pops$pop <- factor( as.character(pops$pop), levels=c("GA", "SC", "WNC", "NC", "KC", "SH", "RB", "OC", "MG", "Mg", "CT", "SR", "BP", "HB", "Hb", "H", "F", "M", "Syc", "NBH", "P", "SLC", "ME", "MDIBL") )
ord <- order(pops$pop)
q <- q[ord,]
plot_data <- q %>%
mutate(id = sample) %>%
gather('cluster', 'prob', V1:V2) %>%
group_by(id) %>%
mutate(likely_assignment = cluster[which.max(prob)],
assingment_prob = max(prob)) %>%
arrange(likely_assignment, desc(assingment_prob)) %>%
ungroup() %>%
dplyr::arrange(factor( as.character(pop), levels=c("GA", "SC", "WNC", "NC", "KC", "SH", "RB", "OC", "MG", "Mg", "CT", "SR", "BP", "HB", "Hb", "H", "F", "M", "Syc", "NBH", "P", "SLC", "ME", "MDIBL") ))
#reorder levels
plot_data$pop <- factor( as.character(plot_data$pop), levels=c("GA", "SC", "WNC", "NC", "KC", "SH", "RB", "OC", "MG", "Mg", "CT", "SR", "BP", "HB", "Hb", "H", "F", "M", "Syc", "NBH", "P", "SLC", "ME", "MDIBL") )
#ggplot(plot_data, aes(id, prob, fill = cluster)) +
# geom_col() +
# theme_classic()
ggplot(plot_data, aes(id, prob, fill = cluster)) +
geom_col() +
facet_grid(~pop, scales = 'free', space = 'free')+ theme(panel.spacing=unit(0.01, "lines"), axis.title.x=element_blank(),axis.text.x=element_blank(),axis.ticks.x=element_blank(), legend.position = "none")
format data for clumpak
(for log in `ls *.log`; do grep -Po 'like=\K[^ ]+' $log; done) > logfile
logs <- as.data.frame(read.table("./ngs_admix/logfile_likelihoods.txt"))
logs$K <- c(rep("1", 10), rep("2",
10), rep("3", 10), rep("4", 10), rep("5", 10), rep("6", 10))
colnames(logs)[1] <- "likelihood"
write.table(logs[, c(2, 1)], "/ngs_admix/logfile_likelihood_formatted.txt", row.names = F,
col.names = F, quote = F)
delta K method suggest best K is 2, but prob(k) is maximixed at k = 6, therefore will run for more K
SFS Analysis Outline
* Make consensus sequence from grandis reads to polarize the SAF - doFasta -2
* estimate SAFs for each population and metapopulation -doSAF * estimate per population SFS
* estimate 2dSFS for each pair of populations
* 1d and 2dSFS used in downstream analyses: per-site summary statistics, global Fst, and dadi
questions: is it possible to estimate SAF without re-estimating GLs (we already have them) what filtering is neccesarry for the population and 2d sfs, when there are so few indiviudals a maf of 1% becomes unmeaningful, should we filter to a higher MAF or at least, increase min number of individuals for a site to be retained (i.e. saf estimate is likely to be very very bad if there are sites with 2 individuals) *how should we pool sampling locales into metapopulations (i.e. it’s clear that all the NBH populations should be combined into a single pop, but what about all NJ populations or all ME populations) - what spatial or genetic scale should be the cutoff for pooling populations
ancestral sequence
Polarize SAF by an outgroup (f grandis) use doFasta -2 with -doCounts 1 and the same filters as $FILTERS from the GL run
#!/bin/bash
# set max wall-clock time (D-HH:MM:SS)
#SBATCH --time=6-23:59:00
#SBATCH --cpus-per-task=38
# set partition/queue to use
#SBATCH --partition=week-long-cpu
# set name of job
#SBATCH --job-name=consensus
# set name of output file
#SBATCH --output=consensus.out
ANGSD="/opt/angsd0.928/angsd/angsd"
REF="/home/ddayan/fundulus/genome/f_het_genome_3.0.2.fa.gz"
FILTERS="-uniqueOnly 1 -remove_bads 1 -trim 0 -C 50 -baq 1 -minMapQ 20"
DOS="-doFasta 2 -doCounts 1 -dumpCounts 2"
$ANGSD -P 38 -b ../meta_data/grandis.bams -ref $REF -out ./Results/grandis_consensus.fa \
$FILTERS $DOS \
SAF 1 big population
require(vegan)
require(ade4)
require(adespatial)
require(SoDA)
require(marmap)
distances
rather than using distance matrices based on latitude and longitude, we make a more biologically relevant model where spatial distances among populations is determined by connections via shallow water <5 meters) this is accomplished via the marmap package in R
atlantic <- getNOAA.bathy(lon1 = -83, lon2 = -65, lat1 = 28, lat2 = 48, resolution = 1)
blues <- c("lightsteelblue4", "lightsteelblue3", "lightsteelblue2", "lightsteelblue1")
greys <- c(grey(0.6), grey(0.93), grey(0.99))
plot(atlantic, image = TRUE, land = TRUE, lwd = 0.03, bpal = list(c(0, max(atlantic), greys),c(min(atlantic), 0, blues)))
#sites<-read.table("~/Science/fundulus_clinal/cline/environmental/lat_lon2.txt", header=T)
sites<-read.table("./environmental/bathy_sites.txt", header = T, sep = "\t")
#remove braytonpt and oc
#sites3<-sites2[c(-7,-11),]
#row.names(sites)<-sites[,1]
#sites<-sites[,-1]
#sites2<-sites
#sites2[,1]<-sites[,2]
#sites2[,2]<-sites[,1]
trans<-trans.mat(atlantic, min.depth=1 , max.depth=-15)
path <- lc.dist(trans, sites[,c(3,2)], res = "path")
## lon lat depth
## 1 -71.18604 41.71250 -5
## 2 -72.52448 41.27902 0
## 3 -81.35452 31.44932 -2
## 4 -71.03655 41.49594 -3
## 5 -76.42462 37.29845 3
## 6 -70.81464 41.65875 -5
## 7 -68.31693 44.44181 41
## 8 -69.72242 43.94622 0
## 9 -74.07045 40.05046 -1
## 10 -70.90760 41.62486 -6
## 11 -75.65578 35.91287 0
## 12 -74.18475 39.80866 1
## 13 -74.32448 39.50878 -3
## 14 -79.91624 32.75873 -5
## 15 -74.75746 39.06786 -3
## 16 -70.97979 41.53583 -1
## 17 -71.52557 41.38235 -6
## 18 -77.92477 34.02835 0
##
|
| | 0%
|
| | 1%
|
|= | 1%
|
|= | 2%
|
|== | 3%
|
|=== | 4%
|
|=== | 5%
|
|==== | 6%
|
|==== | 7%
|
|===== | 7%
|
|===== | 8%
|
|====== | 8%
|
|====== | 9%
|
|====== | 10%
|
|======= | 10%
|
|======= | 11%
|
|======== | 12%
|
|======== | 13%
|
|========= | 14%
|
|========== | 15%
|
|========== | 16%
|
|=========== | 16%
|
|=========== | 17%
|
|=========== | 18%
|
|============ | 18%
|
|============ | 19%
|
|============= | 20%
|
|============== | 21%
|
|============== | 22%
|
|=============== | 23%
|
|=============== | 24%
|
|================ | 24%
|
|================ | 25%
|
|================= | 25%
|
|================= | 26%
|
|================= | 27%
|
|================== | 27%
|
|================== | 28%
|
|=================== | 29%
|
|==================== | 30%
|
|==================== | 31%
|
|===================== | 32%
|
|===================== | 33%
|
|====================== | 33%
|
|====================== | 34%
|
|======================= | 35%
|
|======================= | 36%
|
|======================== | 37%
|
|========================= | 38%
|
|========================= | 39%
|
|========================== | 40%
|
|========================== | 41%
|
|=========================== | 41%
|
|=========================== | 42%
|
|============================ | 42%
|
|============================ | 43%
|
|============================ | 44%
|
|============================= | 44%
|
|============================= | 45%
|
|============================== | 46%
|
|=============================== | 47%
|
|=============================== | 48%
|
|================================ | 49%
|
|================================ | 50%
|
|================================= | 50%
|
|================================= | 51%
|
|================================== | 52%
|
|================================== | 53%
|
|=================================== | 54%
|
|==================================== | 55%
|
|==================================== | 56%
|
|===================================== | 56%
|
|===================================== | 57%
|
|===================================== | 58%
|
|====================================== | 58%
|
|====================================== | 59%
|
|======================================= | 59%
|
|======================================= | 60%
|
|======================================== | 61%
|
|======================================== | 62%
|
|========================================= | 63%
|
|========================================== | 64%
|
|========================================== | 65%
|
|=========================================== | 66%
|
|=========================================== | 67%
|
|============================================ | 67%
|
|============================================ | 68%
|
|============================================= | 69%
|
|============================================= | 70%
|
|============================================== | 71%
|
|=============================================== | 72%
|
|=============================================== | 73%
|
|================================================ | 73%
|
|================================================ | 74%
|
|================================================ | 75%
|
|================================================= | 75%
|
|================================================= | 76%
|
|================================================== | 76%
|
|================================================== | 77%
|
|=================================================== | 78%
|
|=================================================== | 79%
|
|==================================================== | 80%
|
|===================================================== | 81%
|
|===================================================== | 82%
|
|====================================================== | 82%
|
|====================================================== | 83%
|
|====================================================== | 84%
|
|======================================================= | 84%
|
|======================================================= | 85%
|
|======================================================== | 86%
|
|========================================================= | 87%
|
|========================================================= | 88%
|
|========================================================== | 89%
|
|========================================================== | 90%
|
|=========================================================== | 90%
|
|=========================================================== | 91%
|
|=========================================================== | 92%
|
|============================================================ | 92%
|
|============================================================ | 93%
|
|============================================================= | 93%
|
|============================================================= | 94%
|
|============================================================== | 95%
|
|============================================================== | 96%
|
|=============================================================== | 97%
|
|================================================================ | 98%
|
|================================================================ | 99%
|
|=================================================================| 99%
|
|=================================================================| 100%
#distance_matrix <- lc.dist(trans, sites[,c(3,2)], res = "dist")
points(sites[,c(3,2)], pch = 21, col = "blue", bg = col2alpha("blue", .9),cex = 1.2)
lapply(path, lines, col = "orange", lwd = 2, lty = 1)->dummy
#get.depth(atlantic, sites[,c(3,2)] , locator = FALSE)
#distance_matrix
#prep site data for dbmem
#merge all NBH inds to a single pop and convert to all caps
pops <- pops %>%
mutate(pop = as.character(pop)) %>%
mutate(pop = replace(pop, pop == "F", "NBH")) %>%
mutate(pop = replace(pop, pop == "P" , "NBH")) %>%
mutate(pop = replace(pop, pop == "Syc" , "NBH")) %>%
mutate(pop = replace(pop, pop == "H", "NBH")) %>%
mutate(pop = replace(pop, pop == "Mg", "MG")) %>%
mutate(pop = replace(pop, pop == "Hb", "HB"))
ind_sites <- pops %>%
left_join(., sites, by = "pop")
## Warning: Column `pop` joining character vector and factor, coercing into
## character vector
distance_matrix <- lc.dist(trans, ind_sites[,c(4,3)], res = "dist")
## lon lat depth
## 1 -71.18604 41.71250 -5
## 2 -71.18604 41.71250 -5
## 3 -71.18604 41.71250 -5
## 4 -71.18604 41.71250 -5
## 5 -71.18604 41.71250 -5
## 6 -71.18604 41.71250 -5
## 7 -71.18604 41.71250 -5
## 8 -71.18604 41.71250 -5
## 9 -71.18604 41.71250 -5
## 10 -71.18604 41.71250 -5
## 11 -71.18604 41.71250 -5
## 12 -71.18604 41.71250 -5
## 13 -71.18604 41.71250 -5
## 14 -71.18604 41.71250 -5
## 15 -71.18604 41.71250 -5
## 16 -71.18604 41.71250 -5
## 17 -71.18604 41.71250 -5
## 18 -71.18604 41.71250 -5
## 19 -71.18604 41.71250 -5
## 20 -71.18604 41.71250 -5
## 21 -71.18604 41.71250 -5
## 22 -71.18604 41.71250 -5
## 23 -71.18604 41.71250 -5
## 24 -71.18604 41.71250 -5
## 25 -71.18604 41.71250 -5
## 26 -71.18604 41.71250 -5
## 27 -71.18604 41.71250 -5
## 28 -71.18604 41.71250 -5
## 29 -71.18604 41.71250 -5
## 30 -71.18604 41.71250 -5
## 31 -71.18604 41.71250 -5
## 32 -71.18604 41.71250 -5
## 33 -71.18604 41.71250 -5
## 34 -71.18604 41.71250 -5
## 35 -71.18604 41.71250 -5
## 36 -71.18604 41.71250 -5
## 37 -71.18604 41.71250 -5
## 38 -71.18604 41.71250 -5
## 39 -71.18604 41.71250 -5
## 40 -71.18604 41.71250 -5
## 41 -71.18604 41.71250 -5
## 42 -71.18604 41.71250 -5
## 43 -71.18604 41.71250 -5
## 44 -71.18604 41.71250 -5
## 45 -71.18604 41.71250 -5
## 46 -71.18604 41.71250 -5
## 47 -71.18604 41.71250 -5
## 48 -71.18604 41.71250 -5
## 49 -71.18604 41.71250 -5
## 50 -71.18604 41.71250 -5
## 51 -72.52448 41.27902 0
## 52 -72.52448 41.27902 0
## 53 -72.52448 41.27902 0
## 54 -72.52448 41.27902 0
## 55 -72.52448 41.27902 0
## 56 -72.52448 41.27902 0
## 57 -72.52448 41.27902 0
## 58 -72.52448 41.27902 0
## 59 -72.52448 41.27902 0
## 60 -72.52448 41.27902 0
## 61 -72.52448 41.27902 0
## 62 -72.52448 41.27902 0
## 63 -72.52448 41.27902 0
## 64 -72.52448 41.27902 0
## 65 -72.52448 41.27902 0
## 66 -72.52448 41.27902 0
## 67 -72.52448 41.27902 0
## 68 -72.52448 41.27902 0
## 69 -72.52448 41.27902 0
## 70 -72.52448 41.27902 0
## 71 -70.90760 41.62486 -6
## 72 -70.90760 41.62486 -6
## 73 -70.90760 41.62486 -6
## 74 -70.90760 41.62486 -6
## 75 -70.90760 41.62486 -6
## 76 -70.90760 41.62486 -6
## 77 -70.90760 41.62486 -6
## 78 -70.90760 41.62486 -6
## 79 -70.90760 41.62486 -6
## 80 -70.90760 41.62486 -6
## 81 -70.90760 41.62486 -6
## 82 -70.90760 41.62486 -6
## 83 -70.90760 41.62486 -6
## 84 -70.90760 41.62486 -6
## 85 -70.90760 41.62486 -6
## 86 -70.90760 41.62486 -6
## 87 -70.90760 41.62486 -6
## 88 -70.90760 41.62486 -6
## 89 -70.90760 41.62486 -6
## 90 -70.90760 41.62486 -6
## 91 -70.90760 41.62486 -6
## 92 -70.90760 41.62486 -6
## 93 -70.90760 41.62486 -6
## 94 -70.90760 41.62486 -6
## 95 -70.90760 41.62486 -6
## 96 -70.90760 41.62486 -6
## 97 -70.90760 41.62486 -6
## 98 -70.90760 41.62486 -6
## 99 -70.90760 41.62486 -6
## 100 -70.90760 41.62486 -6
## 101 -81.35452 31.44932 -2
## 102 -81.35452 31.44932 -2
## 103 -81.35452 31.44932 -2
## 104 -81.35452 31.44932 -2
## 105 -81.35452 31.44932 -2
## 106 -81.35452 31.44932 -2
## 107 -81.35452 31.44932 -2
## 108 -81.35452 31.44932 -2
## 109 -81.35452 31.44932 -2
## 110 -81.35452 31.44932 -2
## 111 -81.35452 31.44932 -2
## 112 -81.35452 31.44932 -2
## 113 -81.35452 31.44932 -2
## 114 -81.35452 31.44932 -2
## 115 -81.35452 31.44932 -2
## 116 -81.35452 31.44932 -2
## 117 -81.35452 31.44932 -2
## 118 -81.35452 31.44932 -2
## 119 -81.35452 31.44932 -2
## 120 -81.35452 31.44932 -2
## 121 -81.35452 31.44932 -2
## 122 -81.35452 31.44932 -2
## 123 -81.35452 31.44932 -2
## 124 -81.35452 31.44932 -2
## 125 -81.35452 31.44932 -2
## 126 -81.35452 31.44932 -2
## 127 -81.35452 31.44932 -2
## 128 -81.35452 31.44932 -2
## 129 -81.35452 31.44932 -2
## 130 -81.35452 31.44932 -2
## 131 -81.35452 31.44932 -2
## 132 -81.35452 31.44932 -2
## 133 -81.35452 31.44932 -2
## 134 -81.35452 31.44932 -2
## 135 -81.35452 31.44932 -2
## 136 -81.35452 31.44932 -2
## 137 -81.35452 31.44932 -2
## 138 -81.35452 31.44932 -2
## 139 -70.90760 41.62486 -6
## 140 -70.90760 41.62486 -6
## 141 -70.90760 41.62486 -6
## 142 -70.90760 41.62486 -6
## 143 -70.90760 41.62486 -6
## 144 -70.90760 41.62486 -6
## 145 -70.90760 41.62486 -6
## 146 -70.90760 41.62486 -6
## 147 -70.90760 41.62486 -6
## 148 -70.90760 41.62486 -6
## 149 -70.90760 41.62486 -6
## 150 -70.90760 41.62486 -6
## 151 -70.90760 41.62486 -6
## 152 -70.90760 41.62486 -6
## 153 -70.90760 41.62486 -6
## 154 -70.90760 41.62486 -6
## 155 -70.90760 41.62486 -6
## 156 -70.90760 41.62486 -6
## 157 -70.90760 41.62486 -6
## 158 -70.90760 41.62486 -6
## 159 -70.90760 41.62486 -6
## 160 -70.90760 41.62486 -6
## 161 -70.90760 41.62486 -6
## 162 -70.90760 41.62486 -6
## 163 -70.90760 41.62486 -6
## 164 -70.90760 41.62486 -6
## 165 -70.90760 41.62486 -6
## 166 -70.90760 41.62486 -6
## 167 -70.90760 41.62486 -6
## 168 -70.90760 41.62486 -6
## 169 -71.03655 41.49594 -3
## 170 -71.03655 41.49594 -3
## 171 -71.03655 41.49594 -3
## 172 -71.03655 41.49594 -3
## 173 -71.03655 41.49594 -3
## 174 -71.03655 41.49594 -3
## 175 -71.03655 41.49594 -3
## 176 -71.03655 41.49594 -3
## 177 -71.03655 41.49594 -3
## 178 -71.03655 41.49594 -3
## 179 -71.03655 41.49594 -3
## 180 -71.03655 41.49594 -3
## 181 -71.03655 41.49594 -3
## 182 -71.03655 41.49594 -3
## 183 -71.03655 41.49594 -3
## 184 -71.03655 41.49594 -3
## 185 -71.03655 41.49594 -3
## 186 -71.03655 41.49594 -3
## 187 -71.03655 41.49594 -3
## 188 -71.03655 41.49594 -3
## 189 -71.03655 41.49594 -3
## 190 -71.03655 41.49594 -3
## 191 -71.03655 41.49594 -3
## 192 -71.03655 41.49594 -3
## 193 -71.03655 41.49594 -3
## 194 -71.03655 41.49594 -3
## 195 -71.03655 41.49594 -3
## 196 -71.03655 41.49594 -3
## 197 -71.03655 41.49594 -3
## 198 -71.03655 41.49594 -3
## 199 -71.03655 41.49594 -3
## 200 -71.03655 41.49594 -3
## 201 -71.03655 41.49594 -3
## 202 -71.03655 41.49594 -3
## 203 -71.03655 41.49594 -3
## 204 -71.03655 41.49594 -3
## 205 -71.03655 41.49594 -3
## 206 -71.03655 41.49594 -3
## 207 -71.03655 41.49594 -3
## 208 -71.03655 41.49594 -3
## 209 -71.03655 41.49594 -3
## 210 -71.03655 41.49594 -3
## 211 -71.03655 41.49594 -3
## 212 -71.03655 41.49594 -3
## 213 -71.03655 41.49594 -3
## 214 -71.03655 41.49594 -3
## 215 -71.03655 41.49594 -3
## 216 -71.03655 41.49594 -3
## 217 -76.42462 37.29845 3
## 218 -76.42462 37.29845 3
## 219 -76.42462 37.29845 3
## 220 -76.42462 37.29845 3
## 221 -76.42462 37.29845 3
## 222 -76.42462 37.29845 3
## 223 -76.42462 37.29845 3
## 224 -76.42462 37.29845 3
## 225 -76.42462 37.29845 3
## 226 -76.42462 37.29845 3
## 227 -76.42462 37.29845 3
## 228 -76.42462 37.29845 3
## 229 -76.42462 37.29845 3
## 230 -76.42462 37.29845 3
## 231 -76.42462 37.29845 3
## 232 -76.42462 37.29845 3
## 233 -76.42462 37.29845 3
## 234 -76.42462 37.29845 3
## 235 -76.42462 37.29845 3
## 236 -76.42462 37.29845 3
## 237 -76.42462 37.29845 3
## 238 -76.42462 37.29845 3
## 239 -76.42462 37.29845 3
## 240 -76.42462 37.29845 3
## 241 -76.42462 37.29845 3
## 242 -70.81464 41.65875 -5
## 243 -70.81464 41.65875 -5
## 244 -70.81464 41.65875 -5
## 245 -70.81464 41.65875 -5
## 246 -70.81464 41.65875 -5
## 247 -70.81464 41.65875 -5
## 248 -70.81464 41.65875 -5
## 249 -70.81464 41.65875 -5
## 250 -70.81464 41.65875 -5
## 251 -70.81464 41.65875 -5
## 252 -70.81464 41.65875 -5
## 253 -70.81464 41.65875 -5
## 254 -70.81464 41.65875 -5
## 255 -70.81464 41.65875 -5
## 256 -70.81464 41.65875 -5
## 257 -70.81464 41.65875 -5
## 258 -70.81464 41.65875 -5
## 259 -70.81464 41.65875 -5
## 260 -70.81464 41.65875 -5
## 261 -70.81464 41.65875 -5
## 262 -70.81464 41.65875 -5
## 263 -70.81464 41.65875 -5
## 264 -70.81464 41.65875 -5
## 265 -70.81464 41.65875 -5
## 266 -70.81464 41.65875 -5
## 267 -70.81464 41.65875 -5
## 268 -70.81464 41.65875 -5
## 269 -70.81464 41.65875 -5
## 270 -70.81464 41.65875 -5
## 271 -70.81464 41.65875 -5
## 272 -70.81464 41.65875 -5
## 273 -68.31693 44.44181 41
## 274 -68.31693 44.44181 41
## 275 -68.31693 44.44181 41
## 276 -68.31693 44.44181 41
## 277 -68.31693 44.44181 41
## 278 -68.31693 44.44181 41
## 279 -68.31693 44.44181 41
## 280 -68.31693 44.44181 41
## 281 -68.31693 44.44181 41
## 282 -68.31693 44.44181 41
## 283 -68.31693 44.44181 41
## 284 -68.31693 44.44181 41
## 285 -68.31693 44.44181 41
## 286 -68.31693 44.44181 41
## 287 -68.31693 44.44181 41
## 288 -68.31693 44.44181 41
## 289 -68.31693 44.44181 41
## 290 -68.31693 44.44181 41
## 291 -68.31693 44.44181 41
## 292 -68.31693 44.44181 41
## 293 -68.31693 44.44181 41
## 294 -68.31693 44.44181 41
## 295 -68.31693 44.44181 41
## 296 -69.72242 43.94622 0
## 297 -69.72242 43.94622 0
## 298 -69.72242 43.94622 0
## 299 -69.72242 43.94622 0
## 300 -69.72242 43.94622 0
## 301 -69.72242 43.94622 0
## 302 -69.72242 43.94622 0
## 303 -69.72242 43.94622 0
## 304 -69.72242 43.94622 0
## 305 -69.72242 43.94622 0
## 306 -69.72242 43.94622 0
## 307 -69.72242 43.94622 0
## 308 -69.72242 43.94622 0
## 309 -69.72242 43.94622 0
## 310 -69.72242 43.94622 0
## 311 -69.72242 43.94622 0
## 312 -69.72242 43.94622 0
## 313 -69.72242 43.94622 0
## 314 -69.72242 43.94622 0
## 315 -69.72242 43.94622 0
## 316 -69.72242 43.94622 0
## 317 -69.72242 43.94622 0
## 318 -69.72242 43.94622 0
## 319 -69.72242 43.94622 0
## 320 -69.72242 43.94622 0
## 321 -69.72242 43.94622 0
## 322 -69.72242 43.94622 0
## 323 -74.07045 40.05046 -1
## 324 -74.07045 40.05046 -1
## 325 -74.07045 40.05046 -1
## 326 -74.07045 40.05046 -1
## 327 -74.07045 40.05046 -1
## 328 -74.07045 40.05046 -1
## 329 -74.07045 40.05046 -1
## 330 -74.07045 40.05046 -1
## 331 -74.07045 40.05046 -1
## 332 -74.07045 40.05046 -1
## 333 -74.07045 40.05046 -1
## 334 -74.07045 40.05046 -1
## 335 -74.07045 40.05046 -1
## 336 -74.07045 40.05046 -1
## 337 -74.07045 40.05046 -1
## 338 -74.07045 40.05046 -1
## 339 -74.07045 40.05046 -1
## 340 -74.07045 40.05046 -1
## 341 -74.07045 40.05046 -1
## 342 -74.07045 40.05046 -1
## 343 -74.07045 40.05046 -1
## 344 -74.07045 40.05046 -1
## 345 -74.07045 40.05046 -1
## 346 -74.07045 40.05046 -1
## 347 -74.07045 40.05046 -1
## 348 -74.07045 40.05046 -1
## 349 -74.07045 40.05046 -1
## 350 -74.07045 40.05046 -1
## 351 -74.07045 40.05046 -1
## 352 -74.07045 40.05046 -1
## 353 -74.07045 40.05046 -1
## 354 -74.07045 40.05046 -1
## 355 -74.07045 40.05046 -1
## 356 -74.07045 40.05046 -1
## 357 -74.07045 40.05046 -1
## 358 -74.07045 40.05046 -1
## 359 -74.07045 40.05046 -1
## 360 -74.07045 40.05046 -1
## 361 -74.07045 40.05046 -1
## 362 -74.07045 40.05046 -1
## 363 -74.07045 40.05046 -1
## 364 -74.07045 40.05046 -1
## 365 -74.07045 40.05046 -1
## 366 -74.07045 40.05046 -1
## 367 -74.07045 40.05046 -1
## 368 -74.07045 40.05046 -1
## 369 -74.07045 40.05046 -1
## 370 -74.07045 40.05046 -1
## 371 -74.07045 40.05046 -1
## 372 -74.07045 40.05046 -1
## 373 -74.07045 40.05046 -1
## 374 -74.07045 40.05046 -1
## 375 -74.07045 40.05046 -1
## 376 -74.07045 40.05046 -1
## 377 -74.07045 40.05046 -1
## 378 -74.07045 40.05046 -1
## 379 -74.07045 40.05046 -1
## 380 -74.07045 40.05046 -1
## 381 -74.07045 40.05046 -1
## 382 -74.07045 40.05046 -1
## 383 -74.07045 40.05046 -1
## 384 -74.07045 40.05046 -1
## 385 -74.07045 40.05046 -1
## 386 -74.07045 40.05046 -1
## 387 -74.07045 40.05046 -1
## 388 -74.07045 40.05046 -1
## 389 -74.07045 40.05046 -1
## 390 -74.07045 40.05046 -1
## 391 -74.07045 40.05046 -1
## 392 -74.07045 40.05046 -1
## 393 -74.07045 40.05046 -1
## 394 -74.07045 40.05046 -1
## 395 -74.07045 40.05046 -1
## 396 -74.07045 40.05046 -1
## 397 -74.07045 40.05046 -1
## 398 -74.07045 40.05046 -1
## 399 -74.07045 40.05046 -1
## 400 -74.07045 40.05046 -1
## 401 -74.07045 40.05046 -1
## 402 -74.07045 40.05046 -1
## 403 -74.07045 40.05046 -1
## 404 -74.07045 40.05046 -1
## 405 -74.07045 40.05046 -1
## 406 -74.07045 40.05046 -1
## 407 -74.07045 40.05046 -1
## 408 -74.07045 40.05046 -1
## 409 -74.07045 40.05046 -1
## 410 -74.07045 40.05046 -1
## 411 -74.07045 40.05046 -1
## 412 -74.07045 40.05046 -1
## 413 -74.07045 40.05046 -1
## 414 -74.07045 40.05046 -1
## 415 -74.07045 40.05046 -1
## 416 -74.07045 40.05046 -1
## 417 -74.07045 40.05046 -1
## 418 -74.07045 40.05046 -1
## 419 -74.07045 40.05046 -1
## 420 -74.07045 40.05046 -1
## 421 -74.07045 40.05046 -1
## 422 -74.07045 40.05046 -1
## 423 -74.07045 40.05046 -1
## 424 -74.07045 40.05046 -1
## 425 -74.07045 40.05046 -1
## 426 -74.07045 40.05046 -1
## 427 -74.07045 40.05046 -1
## 428 -74.07045 40.05046 -1
## 429 -74.07045 40.05046 -1
## 430 -74.07045 40.05046 -1
## 431 -74.07045 40.05046 -1
## 432 -74.07045 40.05046 -1
## 433 -74.07045 40.05046 -1
## 434 -74.07045 40.05046 -1
## 435 -74.07045 40.05046 -1
## 436 -74.07045 40.05046 -1
## 437 -74.07045 40.05046 -1
## 438 -74.07045 40.05046 -1
## 439 -74.07045 40.05046 -1
## 440 -74.07045 40.05046 -1
## 441 -74.07045 40.05046 -1
## 442 -74.07045 40.05046 -1
## 443 -74.07045 40.05046 -1
## 444 -74.07045 40.05046 -1
## 445 -74.07045 40.05046 -1
## 446 -74.07045 40.05046 -1
## 447 -74.07045 40.05046 -1
## 448 -74.07045 40.05046 -1
## 449 -74.07045 40.05046 -1
## 450 -74.07045 40.05046 -1
## 451 -74.07045 40.05046 -1
## 452 -74.07045 40.05046 -1
## 453 -74.07045 40.05046 -1
## 454 -74.07045 40.05046 -1
## 455 -74.07045 40.05046 -1
## 456 -74.07045 40.05046 -1
## 457 -74.07045 40.05046 -1
## 458 -74.07045 40.05046 -1
## 459 -74.07045 40.05046 -1
## 460 -74.07045 40.05046 -1
## 461 -74.07045 40.05046 -1
## 462 -74.07045 40.05046 -1
## 463 -74.07045 40.05046 -1
## 464 -74.07045 40.05046 -1
## 465 -74.07045 40.05046 -1
## 466 -74.07045 40.05046 -1
## 467 -74.07045 40.05046 -1
## 468 -74.07045 40.05046 -1
## 469 -74.07045 40.05046 -1
## 470 -74.07045 40.05046 -1
## 471 -74.07045 40.05046 -1
## 472 -74.07045 40.05046 -1
## 473 -74.07045 40.05046 -1
## 474 -74.07045 40.05046 -1
## 475 -74.07045 40.05046 -1
## 476 -74.07045 40.05046 -1
## 477 -74.07045 40.05046 -1
## 478 -74.07045 40.05046 -1
## 479 -74.07045 40.05046 -1
## 480 -74.07045 40.05046 -1
## 481 -74.07045 40.05046 -1
## 482 -74.07045 40.05046 -1
## 483 -74.07045 40.05046 -1
## 484 -74.07045 40.05046 -1
## 485 -74.07045 40.05046 -1
## 486 -74.07045 40.05046 -1
## 487 -74.07045 40.05046 -1
## 488 -74.07045 40.05046 -1
## 489 -74.07045 40.05046 -1
## 490 -74.07045 40.05046 -1
## 491 -74.07045 40.05046 -1
## 492 -74.07045 40.05046 -1
## 493 -74.07045 40.05046 -1
## 494 -74.07045 40.05046 -1
## 495 -74.07045 40.05046 -1
## 496 -74.07045 40.05046 -1
## 497 -74.07045 40.05046 -1
## 498 -74.07045 40.05046 -1
## 499 -74.07045 40.05046 -1
## 500 -74.07045 40.05046 -1
## 501 -74.07045 40.05046 -1
## 502 -74.07045 40.05046 -1
## 503 -74.07045 40.05046 -1
## 504 -74.07045 40.05046 -1
## 505 -74.07045 40.05046 -1
## 506 -74.07045 40.05046 -1
## 507 -74.07045 40.05046 -1
## 508 -74.07045 40.05046 -1
## 509 -74.07045 40.05046 -1
## 510 -74.07045 40.05046 -1
## 511 -74.07045 40.05046 -1
## 512 -74.07045 40.05046 -1
## 513 -74.07045 40.05046 -1
## 514 -74.07045 40.05046 -1
## 515 -74.07045 40.05046 -1
## 516 -74.07045 40.05046 -1
## 517 -74.07045 40.05046 -1
## 518 -74.07045 40.05046 -1
## 519 -74.07045 40.05046 -1
## 520 -74.07045 40.05046 -1
## 521 -74.07045 40.05046 -1
## 522 -74.07045 40.05046 -1
## 523 -74.07045 40.05046 -1
## 524 -74.07045 40.05046 -1
## 525 -74.07045 40.05046 -1
## 526 -74.07045 40.05046 -1
## 527 -74.07045 40.05046 -1
## 528 -74.07045 40.05046 -1
## 529 -74.07045 40.05046 -1
## 530 -74.07045 40.05046 -1
## 531 -74.07045 40.05046 -1
## 532 -74.07045 40.05046 -1
## 533 -74.07045 40.05046 -1
## 534 -74.07045 40.05046 -1
## 535 -74.07045 40.05046 -1
## 536 -74.07045 40.05046 -1
## 537 -74.07045 40.05046 -1
## 538 -74.07045 40.05046 -1
## 539 -74.07045 40.05046 -1
## 540 -74.07045 40.05046 -1
## 541 -74.07045 40.05046 -1
## 542 -74.07045 40.05046 -1
## 543 -74.07045 40.05046 -1
## 544 -74.07045 40.05046 -1
## 545 -74.07045 40.05046 -1
## 546 -74.07045 40.05046 -1
## 547 -74.07045 40.05046 -1
## 548 -74.07045 40.05046 -1
## 549 -74.07045 40.05046 -1
## 550 -74.07045 40.05046 -1
## 551 -74.07045 40.05046 -1
## 552 -74.07045 40.05046 -1
## 553 -74.07045 40.05046 -1
## 554 -74.07045 40.05046 -1
## 555 -74.07045 40.05046 -1
## 556 -74.07045 40.05046 -1
## 557 -74.07045 40.05046 -1
## 558 -74.07045 40.05046 -1
## 559 -74.07045 40.05046 -1
## 560 -74.07045 40.05046 -1
## 561 -74.07045 40.05046 -1
## 562 -74.07045 40.05046 -1
## 563 -74.07045 40.05046 -1
## 564 -74.07045 40.05046 -1
## 565 -74.07045 40.05046 -1
## 566 -74.07045 40.05046 -1
## 567 -74.07045 40.05046 -1
## 568 -74.07045 40.05046 -1
## 569 -74.07045 40.05046 -1
## 570 -74.07045 40.05046 -1
## 571 -74.07045 40.05046 -1
## 572 -74.07045 40.05046 -1
## 573 -74.07045 40.05046 -1
## 574 -74.07045 40.05046 -1
## 575 -74.07045 40.05046 -1
## 576 -74.07045 40.05046 -1
## 577 -74.07045 40.05046 -1
## 578 -74.07045 40.05046 -1
## 579 -74.07045 40.05046 -1
## 580 -74.07045 40.05046 -1
## 581 -74.07045 40.05046 -1
## 582 -74.07045 40.05046 -1
## 583 -74.07045 40.05046 -1
## 584 -74.07045 40.05046 -1
## 585 -74.07045 40.05046 -1
## 586 -74.07045 40.05046 -1
## 587 -74.07045 40.05046 -1
## 588 -74.07045 40.05046 -1
## 589 -74.07045 40.05046 -1
## 590 -74.07045 40.05046 -1
## 591 -74.07045 40.05046 -1
## 592 -74.07045 40.05046 -1
## 593 -74.07045 40.05046 -1
## 594 -74.07045 40.05046 -1
## 595 -74.07045 40.05046 -1
## 596 -74.07045 40.05046 -1
## 597 -74.07045 40.05046 -1
## 598 -74.07045 40.05046 -1
## 599 -74.07045 40.05046 -1
## 600 -74.07045 40.05046 -1
## 601 -74.07045 40.05046 -1
## 602 -74.07045 40.05046 -1
## 603 -74.07045 40.05046 -1
## 604 -74.07045 40.05046 -1
## 605 -74.07045 40.05046 -1
## 606 -74.07045 40.05046 -1
## 607 -74.07045 40.05046 -1
## 608 -74.07045 40.05046 -1
## 609 -74.07045 40.05046 -1
## 610 -74.07045 40.05046 -1
## 611 -74.07045 40.05046 -1
## 612 -74.07045 40.05046 -1
## 613 -74.07045 40.05046 -1
## 614 -74.07045 40.05046 -1
## 615 -74.07045 40.05046 -1
## 616 -74.07045 40.05046 -1
## 617 -74.07045 40.05046 -1
## 618 -74.07045 40.05046 -1
## 619 -74.07045 40.05046 -1
## 620 -74.07045 40.05046 -1
## 621 -74.07045 40.05046 -1
## 622 -74.07045 40.05046 -1
## 623 -74.07045 40.05046 -1
## 624 -74.07045 40.05046 -1
## 625 -74.07045 40.05046 -1
## 626 -74.07045 40.05046 -1
## 627 -74.07045 40.05046 -1
## 628 -74.07045 40.05046 -1
## 629 -74.07045 40.05046 -1
## 630 -74.07045 40.05046 -1
## 631 -74.07045 40.05046 -1
## 632 -74.07045 40.05046 -1
## 633 -74.07045 40.05046 -1
## 634 -74.07045 40.05046 -1
## 635 -74.07045 40.05046 -1
## 636 -74.07045 40.05046 -1
## 637 -70.90760 41.62486 -6
## 638 -70.90760 41.62486 -6
## 639 -70.90760 41.62486 -6
## 640 -70.90760 41.62486 -6
## 641 -70.90760 41.62486 -6
## 642 -70.90760 41.62486 -6
## 643 -70.90760 41.62486 -6
## 644 -70.90760 41.62486 -6
## 645 -70.90760 41.62486 -6
## 646 -70.90760 41.62486 -6
## 647 -70.90760 41.62486 -6
## 648 -70.90760 41.62486 -6
## 649 -70.90760 41.62486 -6
## 650 -70.90760 41.62486 -6
## 651 -70.90760 41.62486 -6
## 652 -70.90760 41.62486 -6
## 653 -70.90760 41.62486 -6
## 654 -70.90760 41.62486 -6
## 655 -70.90760 41.62486 -6
## 656 -70.90760 41.62486 -6
## 657 -70.90760 41.62486 -6
## 658 -70.90760 41.62486 -6
## 659 -70.90760 41.62486 -6
## 660 -70.90760 41.62486 -6
## 661 -70.90760 41.62486 -6
## 662 -70.90760 41.62486 -6
## 663 -70.90760 41.62486 -6
## 664 -70.90760 41.62486 -6
## 665 -70.90760 41.62486 -6
## 666 -70.90760 41.62486 -6
## 667 -75.65578 35.91287 0
## 668 -75.65578 35.91287 0
## 669 -75.65578 35.91287 0
## 670 -75.65578 35.91287 0
## 671 -75.65578 35.91287 0
## 672 -75.65578 35.91287 0
## 673 -75.65578 35.91287 0
## 674 -75.65578 35.91287 0
## 675 -75.65578 35.91287 0
## 676 -75.65578 35.91287 0
## 677 -75.65578 35.91287 0
## 678 -75.65578 35.91287 0
## 679 -75.65578 35.91287 0
## 680 -75.65578 35.91287 0
## 681 -75.65578 35.91287 0
## 682 -75.65578 35.91287 0
## 683 -75.65578 35.91287 0
## 684 -75.65578 35.91287 0
## 685 -75.65578 35.91287 0
## 686 -75.65578 35.91287 0
## 687 -75.65578 35.91287 0
## 688 -75.65578 35.91287 0
## 689 -74.18475 39.80866 1
## 690 -74.18475 39.80866 1
## 691 -74.18475 39.80866 1
## 692 -74.18475 39.80866 1
## 693 -74.18475 39.80866 1
## 694 -74.18475 39.80866 1
## 695 -74.18475 39.80866 1
## 696 -74.18475 39.80866 1
## 697 -74.18475 39.80866 1
## 698 -74.18475 39.80866 1
## 699 -74.18475 39.80866 1
## 700 -74.18475 39.80866 1
## 701 -74.18475 39.80866 1
## 702 -74.18475 39.80866 1
## 703 -74.18475 39.80866 1
## 704 -74.18475 39.80866 1
## 705 -74.18475 39.80866 1
## 706 -74.18475 39.80866 1
## 707 -74.18475 39.80866 1
## 708 -74.18475 39.80866 1
## 709 -74.18475 39.80866 1
## 710 -74.18475 39.80866 1
## 711 -74.18475 39.80866 1
## 712 -74.18475 39.80866 1
## 713 -74.18475 39.80866 1
## 714 -74.18475 39.80866 1
## 715 -74.18475 39.80866 1
## 716 -74.18475 39.80866 1
## 717 -74.18475 39.80866 1
## 718 -74.18475 39.80866 1
## 719 -74.18475 39.80866 1
## 720 -74.18475 39.80866 1
## 721 -74.18475 39.80866 1
## 722 -74.18475 39.80866 1
## 723 -74.18475 39.80866 1
## 724 -74.18475 39.80866 1
## 725 -74.18475 39.80866 1
## 726 -74.18475 39.80866 1
## 727 -74.18475 39.80866 1
## 728 -74.18475 39.80866 1
## 729 -74.18475 39.80866 1
## 730 -74.18475 39.80866 1
## 731 -74.18475 39.80866 1
## 732 -70.90760 41.62486 -6
## 733 -70.90760 41.62486 -6
## 734 -70.90760 41.62486 -6
## 735 -70.90760 41.62486 -6
## 736 -70.90760 41.62486 -6
## 737 -70.90760 41.62486 -6
## 738 -70.90760 41.62486 -6
## 739 -70.90760 41.62486 -6
## 740 -70.90760 41.62486 -6
## 741 -70.90760 41.62486 -6
## 742 -70.90760 41.62486 -6
## 743 -70.90760 41.62486 -6
## 744 -70.90760 41.62486 -6
## 745 -70.90760 41.62486 -6
## 746 -70.90760 41.62486 -6
## 747 -70.90760 41.62486 -6
## 748 -70.90760 41.62486 -6
## 749 -70.90760 41.62486 -6
## 750 -70.90760 41.62486 -6
## 751 -70.90760 41.62486 -6
## 752 -70.90760 41.62486 -6
## 753 -70.90760 41.62486 -6
## 754 -70.90760 41.62486 -6
## 755 -70.90760 41.62486 -6
## 756 -70.90760 41.62486 -6
## 757 -70.90760 41.62486 -6
## 758 -70.90760 41.62486 -6
## 759 -70.90760 41.62486 -6
## 760 -70.90760 41.62486 -6
## 761 -70.90760 41.62486 -6
## 762 -74.32448 39.50878 -3
## 763 -74.32448 39.50878 -3
## 764 -74.32448 39.50878 -3
## 765 -74.32448 39.50878 -3
## 766 -74.32448 39.50878 -3
## 767 -74.32448 39.50878 -3
## 768 -74.32448 39.50878 -3
## 769 -74.32448 39.50878 -3
## 770 -74.32448 39.50878 -3
## 771 -74.32448 39.50878 -3
## 772 -74.32448 39.50878 -3
## 773 -74.32448 39.50878 -3
## 774 -74.32448 39.50878 -3
## 775 -74.32448 39.50878 -3
## 776 -74.32448 39.50878 -3
## 777 -74.32448 39.50878 -3
## 778 -74.32448 39.50878 -3
## 779 -74.32448 39.50878 -3
## 780 -74.32448 39.50878 -3
## 781 -74.32448 39.50878 -3
## 782 -74.32448 39.50878 -3
## 783 -74.32448 39.50878 -3
## 784 -74.32448 39.50878 -3
## 785 -74.32448 39.50878 -3
## 786 -74.32448 39.50878 -3
## 787 -74.32448 39.50878 -3
## 788 -74.32448 39.50878 -3
## 789 -74.32448 39.50878 -3
## 790 -74.32448 39.50878 -3
## 791 -74.32448 39.50878 -3
## 792 -74.32448 39.50878 -3
## 793 -74.32448 39.50878 -3
## 794 -74.32448 39.50878 -3
## 795 -74.32448 39.50878 -3
## 796 -74.32448 39.50878 -3
## 797 -74.32448 39.50878 -3
## 798 -74.32448 39.50878 -3
## 799 -74.32448 39.50878 -3
## 800 -74.32448 39.50878 -3
## 801 -74.32448 39.50878 -3
## 802 -74.32448 39.50878 -3
## 803 -74.32448 39.50878 -3
## 804 -74.32448 39.50878 -3
## 805 -74.32448 39.50878 -3
## 806 -74.32448 39.50878 -3
## 807 -74.32448 39.50878 -3
## 808 -74.32448 39.50878 -3
## 809 -74.32448 39.50878 -3
## 810 -74.32448 39.50878 -3
## 811 -74.32448 39.50878 -3
## 812 -74.32448 39.50878 -3
## 813 -74.32448 39.50878 -3
## 814 -74.32448 39.50878 -3
## 815 -74.32448 39.50878 -3
## 816 -74.32448 39.50878 -3
## 817 -74.32448 39.50878 -3
## 818 -74.32448 39.50878 -3
## 819 -74.32448 39.50878 -3
## 820 -74.32448 39.50878 -3
## 821 -74.32448 39.50878 -3
## 822 -74.32448 39.50878 -3
## 823 -74.32448 39.50878 -3
## 824 -74.32448 39.50878 -3
## 825 -74.32448 39.50878 -3
## 826 -74.32448 39.50878 -3
## 827 -74.32448 39.50878 -3
## 828 -74.32448 39.50878 -3
## 829 -74.32448 39.50878 -3
## 830 -74.32448 39.50878 -3
## 831 -74.32448 39.50878 -3
## 832 -74.32448 39.50878 -3
## 833 -74.32448 39.50878 -3
## 834 -74.32448 39.50878 -3
## 835 -74.32448 39.50878 -3
## 836 -74.32448 39.50878 -3
## 837 -74.32448 39.50878 -3
## 838 -74.32448 39.50878 -3
## 839 -74.32448 39.50878 -3
## 840 -74.32448 39.50878 -3
## 841 -74.32448 39.50878 -3
## 842 -74.32448 39.50878 -3
## 843 -74.32448 39.50878 -3
## 844 -74.32448 39.50878 -3
## 845 -74.32448 39.50878 -3
## 846 -74.32448 39.50878 -3
## 847 -74.32448 39.50878 -3
## 848 -74.32448 39.50878 -3
## 849 -74.32448 39.50878 -3
## 850 -74.32448 39.50878 -3
## 851 -74.32448 39.50878 -3
## 852 -74.32448 39.50878 -3
## 853 -74.32448 39.50878 -3
## 854 -74.32448 39.50878 -3
## 855 -74.32448 39.50878 -3
## 856 -74.32448 39.50878 -3
## 857 -74.32448 39.50878 -3
## 858 -74.32448 39.50878 -3
## 859 -74.32448 39.50878 -3
## 860 -74.32448 39.50878 -3
## 861 -74.32448 39.50878 -3
## 862 -74.32448 39.50878 -3
## 863 -74.32448 39.50878 -3
## 864 -74.32448 39.50878 -3
## 865 -74.32448 39.50878 -3
## 866 -74.32448 39.50878 -3
## 867 -74.32448 39.50878 -3
## 868 -74.32448 39.50878 -3
## 869 -74.32448 39.50878 -3
## 870 -74.32448 39.50878 -3
## 871 -74.32448 39.50878 -3
## 872 -74.32448 39.50878 -3
## 873 -74.32448 39.50878 -3
## 874 -74.32448 39.50878 -3
## 875 -74.32448 39.50878 -3
## 876 -74.32448 39.50878 -3
## 877 -74.32448 39.50878 -3
## 878 -74.32448 39.50878 -3
## 879 -74.32448 39.50878 -3
## 880 -74.32448 39.50878 -3
## 881 -74.32448 39.50878 -3
## 882 -74.32448 39.50878 -3
## 883 -74.32448 39.50878 -3
## 884 -74.32448 39.50878 -3
## 885 -74.32448 39.50878 -3
## 886 -74.32448 39.50878 -3
## 887 -74.32448 39.50878 -3
## 888 -74.32448 39.50878 -3
## 889 -74.32448 39.50878 -3
## 890 -74.32448 39.50878 -3
## 891 -74.32448 39.50878 -3
## 892 -74.32448 39.50878 -3
## 893 -74.32448 39.50878 -3
## 894 -74.32448 39.50878 -3
## 895 -74.32448 39.50878 -3
## 896 -74.32448 39.50878 -3
## 897 -74.32448 39.50878 -3
## 898 -74.32448 39.50878 -3
## 899 -74.32448 39.50878 -3
## 900 -74.32448 39.50878 -3
## 901 -74.32448 39.50878 -3
## 902 -74.32448 39.50878 -3
## 903 -74.32448 39.50878 -3
## 904 -74.32448 39.50878 -3
## 905 -74.32448 39.50878 -3
## 906 -74.32448 39.50878 -3
## 907 -74.32448 39.50878 -3
## 908 -74.32448 39.50878 -3
## 909 -74.32448 39.50878 -3
## 910 -74.32448 39.50878 -3
## 911 -74.32448 39.50878 -3
## 912 -74.32448 39.50878 -3
## 913 -74.32448 39.50878 -3
## 914 -74.32448 39.50878 -3
## 915 -74.32448 39.50878 -3
## 916 -74.32448 39.50878 -3
## 917 -74.32448 39.50878 -3
## 918 -74.32448 39.50878 -3
## 919 -74.32448 39.50878 -3
## 920 -74.32448 39.50878 -3
## 921 -74.32448 39.50878 -3
## 922 -74.32448 39.50878 -3
## 923 -74.32448 39.50878 -3
## 924 -74.32448 39.50878 -3
## 925 -74.32448 39.50878 -3
## 926 -74.32448 39.50878 -3
## 927 -74.32448 39.50878 -3
## 928 -74.32448 39.50878 -3
## 929 -74.32448 39.50878 -3
## 930 -74.32448 39.50878 -3
## 931 -74.32448 39.50878 -3
## 932 -74.32448 39.50878 -3
## 933 -74.32448 39.50878 -3
## 934 -74.32448 39.50878 -3
## 935 -74.32448 39.50878 -3
## 936 -74.32448 39.50878 -3
## 937 -74.32448 39.50878 -3
## 938 -74.32448 39.50878 -3
## 939 -74.32448 39.50878 -3
## 940 -74.32448 39.50878 -3
## 941 -74.32448 39.50878 -3
## 942 -74.32448 39.50878 -3
## 943 -74.32448 39.50878 -3
## 944 -74.32448 39.50878 -3
## 945 -74.32448 39.50878 -3
## 946 -74.32448 39.50878 -3
## 947 -74.32448 39.50878 -3
## 948 -74.32448 39.50878 -3
## 949 -74.32448 39.50878 -3
## 950 -74.32448 39.50878 -3
## 951 -74.32448 39.50878 -3
## 952 -74.32448 39.50878 -3
## 953 -74.32448 39.50878 -3
## 954 -74.32448 39.50878 -3
## 955 -74.32448 39.50878 -3
## 956 -74.32448 39.50878 -3
## 957 -74.32448 39.50878 -3
## 958 -74.32448 39.50878 -3
## 959 -74.32448 39.50878 -3
## 960 -74.32448 39.50878 -3
## 961 -74.32448 39.50878 -3
## 962 -74.32448 39.50878 -3
## 963 -74.32448 39.50878 -3
## 964 -74.32448 39.50878 -3
## 965 -74.32448 39.50878 -3
## 966 -74.32448 39.50878 -3
## 967 -74.32448 39.50878 -3
## 968 -74.32448 39.50878 -3
## 969 -74.32448 39.50878 -3
## 970 -74.32448 39.50878 -3
## 971 -74.32448 39.50878 -3
## 972 -74.32448 39.50878 -3
## 973 -74.32448 39.50878 -3
## 974 -74.32448 39.50878 -3
## 975 -74.32448 39.50878 -3
## 976 -74.32448 39.50878 -3
## 977 -74.32448 39.50878 -3
## 978 -74.32448 39.50878 -3
## 979 -74.32448 39.50878 -3
## 980 -74.32448 39.50878 -3
## 981 -74.32448 39.50878 -3
## 982 -74.32448 39.50878 -3
## 983 -74.32448 39.50878 -3
## 984 -74.32448 39.50878 -3
## 985 -74.32448 39.50878 -3
## 986 -74.32448 39.50878 -3
## 987 -74.32448 39.50878 -3
## 988 -74.32448 39.50878 -3
## 989 -74.32448 39.50878 -3
## 990 -74.32448 39.50878 -3
## 991 -74.32448 39.50878 -3
## 992 -74.32448 39.50878 -3
## 993 -74.32448 39.50878 -3
## 994 -74.32448 39.50878 -3
## 995 -74.32448 39.50878 -3
## 996 -74.32448 39.50878 -3
## 997 -74.32448 39.50878 -3
## 998 -74.32448 39.50878 -3
## 999 -74.32448 39.50878 -3
## 1000 -74.32448 39.50878 -3
## 1001 -74.32448 39.50878 -3
## 1002 -74.32448 39.50878 -3
## 1003 -74.32448 39.50878 -3
## 1004 -74.32448 39.50878 -3
## 1005 -74.32448 39.50878 -3
## 1006 -74.32448 39.50878 -3
## 1007 -74.32448 39.50878 -3
## 1008 -74.32448 39.50878 -3
## 1009 -74.32448 39.50878 -3
## 1010 -74.32448 39.50878 -3
## 1011 -74.32448 39.50878 -3
## 1012 -74.32448 39.50878 -3
## 1013 -74.32448 39.50878 -3
## 1014 -74.32448 39.50878 -3
## 1015 -74.32448 39.50878 -3
## 1016 -74.32448 39.50878 -3
## 1017 -74.32448 39.50878 -3
## 1018 -74.32448 39.50878 -3
## 1019 -74.32448 39.50878 -3
## 1020 -74.32448 39.50878 -3
## 1021 -74.32448 39.50878 -3
## 1022 -74.32448 39.50878 -3
## 1023 -74.32448 39.50878 -3
## 1024 -74.32448 39.50878 -3
## 1025 -74.32448 39.50878 -3
## 1026 -74.32448 39.50878 -3
## 1027 -74.32448 39.50878 -3
## 1028 -74.32448 39.50878 -3
## 1029 -74.32448 39.50878 -3
## 1030 -74.32448 39.50878 -3
## 1031 -74.32448 39.50878 -3
## 1032 -74.32448 39.50878 -3
## 1033 -74.32448 39.50878 -3
## 1034 -74.32448 39.50878 -3
## 1035 -74.32448 39.50878 -3
## 1036 -74.32448 39.50878 -3
## 1037 -74.32448 39.50878 -3
## 1038 -74.32448 39.50878 -3
## 1039 -74.32448 39.50878 -3
## 1040 -74.32448 39.50878 -3
## 1041 -74.32448 39.50878 -3
## 1042 -74.32448 39.50878 -3
## 1043 -74.32448 39.50878 -3
## 1044 -74.32448 39.50878 -3
## 1045 -74.32448 39.50878 -3
## 1046 -74.32448 39.50878 -3
## 1047 -74.32448 39.50878 -3
## 1048 -74.32448 39.50878 -3
## 1049 -74.32448 39.50878 -3
## 1050 -74.32448 39.50878 -3
## 1051 -74.32448 39.50878 -3
## 1052 -74.32448 39.50878 -3
## 1053 -74.32448 39.50878 -3
## 1054 -74.32448 39.50878 -3
## 1055 -74.32448 39.50878 -3
## 1056 -74.32448 39.50878 -3
## 1057 -74.32448 39.50878 -3
## 1058 -74.32448 39.50878 -3
## 1059 -74.32448 39.50878 -3
## 1060 -74.32448 39.50878 -3
## 1061 -74.32448 39.50878 -3
## 1062 -74.32448 39.50878 -3
## 1063 -74.32448 39.50878 -3
## 1064 -74.32448 39.50878 -3
## 1065 -74.32448 39.50878 -3
## 1066 -74.32448 39.50878 -3
## 1067 -74.32448 39.50878 -3
## 1068 -74.32448 39.50878 -3
## 1069 -74.32448 39.50878 -3
## 1070 -74.32448 39.50878 -3
## 1071 -74.32448 39.50878 -3
## 1072 -74.32448 39.50878 -3
## 1073 -74.32448 39.50878 -3
## 1074 -74.32448 39.50878 -3
## 1075 -74.32448 39.50878 -3
## 1076 -74.32448 39.50878 -3
## 1077 -74.32448 39.50878 -3
## 1078 -74.32448 39.50878 -3
## 1079 -74.32448 39.50878 -3
## 1080 -74.32448 39.50878 -3
## 1081 -74.32448 39.50878 -3
## 1082 -74.32448 39.50878 -3
## 1083 -74.32448 39.50878 -3
## 1084 -74.32448 39.50878 -3
## 1085 -74.32448 39.50878 -3
## 1086 -74.32448 39.50878 -3
## 1087 -74.32448 39.50878 -3
## 1088 -74.32448 39.50878 -3
## 1089 -74.32448 39.50878 -3
## 1090 -74.32448 39.50878 -3
## 1091 -74.32448 39.50878 -3
## 1092 -74.32448 39.50878 -3
## 1093 -74.32448 39.50878 -3
## 1094 -74.32448 39.50878 -3
## 1095 -74.32448 39.50878 -3
## 1096 -74.32448 39.50878 -3
## 1097 -74.32448 39.50878 -3
## 1098 -74.32448 39.50878 -3
## 1099 -74.32448 39.50878 -3
## 1100 -74.32448 39.50878 -3
## 1101 -74.32448 39.50878 -3
## 1102 -74.32448 39.50878 -3
## 1103 -74.32448 39.50878 -3
## 1104 -74.32448 39.50878 -3
## 1105 -74.32448 39.50878 -3
## 1106 -74.32448 39.50878 -3
## 1107 -74.32448 39.50878 -3
## 1108 -74.32448 39.50878 -3
## 1109 -74.32448 39.50878 -3
## 1110 -74.32448 39.50878 -3
## 1111 -74.32448 39.50878 -3
## 1112 -74.32448 39.50878 -3
## 1113 -74.32448 39.50878 -3
## 1114 -74.32448 39.50878 -3
## 1115 -74.32448 39.50878 -3
## 1116 -74.32448 39.50878 -3
## 1117 -74.32448 39.50878 -3
## 1118 -74.32448 39.50878 -3
## 1119 -74.32448 39.50878 -3
## 1120 -74.32448 39.50878 -3
## 1121 -74.32448 39.50878 -3
## 1122 -74.32448 39.50878 -3
## 1123 -74.32448 39.50878 -3
## 1124 -74.32448 39.50878 -3
## 1125 -74.32448 39.50878 -3
## 1126 -74.32448 39.50878 -3
## 1127 -74.32448 39.50878 -3
## 1128 -74.32448 39.50878 -3
## 1129 -74.32448 39.50878 -3
## 1130 -74.32448 39.50878 -3
## 1131 -74.32448 39.50878 -3
## 1132 -74.32448 39.50878 -3
## 1133 -74.32448 39.50878 -3
## 1134 -74.32448 39.50878 -3
## 1135 -74.32448 39.50878 -3
## 1136 -74.32448 39.50878 -3
## 1137 -74.32448 39.50878 -3
## 1138 -74.32448 39.50878 -3
## 1139 -74.32448 39.50878 -3
## 1140 -74.32448 39.50878 -3
## 1141 -74.32448 39.50878 -3
## 1142 -74.32448 39.50878 -3
## 1143 -74.32448 39.50878 -3
## 1144 -74.32448 39.50878 -3
## 1145 -74.32448 39.50878 -3
## 1146 -74.32448 39.50878 -3
## 1147 -74.32448 39.50878 -3
## 1148 -74.32448 39.50878 -3
## 1149 -74.32448 39.50878 -3
## 1150 -74.32448 39.50878 -3
## 1151 -74.32448 39.50878 -3
## 1152 -74.32448 39.50878 -3
## 1153 -74.32448 39.50878 -3
## 1154 -79.91624 32.75873 -5
## 1155 -79.91624 32.75873 -5
## 1156 -79.91624 32.75873 -5
## 1157 -79.91624 32.75873 -5
## 1158 -79.91624 32.75873 -5
## 1159 -79.91624 32.75873 -5
## 1160 -79.91624 32.75873 -5
## 1161 -79.91624 32.75873 -5
## 1162 -79.91624 32.75873 -5
## 1163 -79.91624 32.75873 -5
## 1164 -79.91624 32.75873 -5
## 1165 -79.91624 32.75873 -5
## 1166 -79.91624 32.75873 -5
## 1167 -79.91624 32.75873 -5
## 1168 -79.91624 32.75873 -5
## 1169 -79.91624 32.75873 -5
## 1170 -79.91624 32.75873 -5
## 1171 -79.91624 32.75873 -5
## 1172 -79.91624 32.75873 -5
## 1173 -79.91624 32.75873 -5
## 1174 -79.91624 32.75873 -5
## 1175 -79.91624 32.75873 -5
## 1176 -79.91624 32.75873 -5
## 1177 -74.75746 39.06786 -3
## 1178 -74.75746 39.06786 -3
## 1179 -74.75746 39.06786 -3
## 1180 -74.75746 39.06786 -3
## 1181 -74.75746 39.06786 -3
## 1182 -74.75746 39.06786 -3
## 1183 -74.75746 39.06786 -3
## 1184 -74.75746 39.06786 -3
## 1185 -74.75746 39.06786 -3
## 1186 -74.75746 39.06786 -3
## 1187 -74.75746 39.06786 -3
## 1188 -74.75746 39.06786 -3
## 1189 -74.75746 39.06786 -3
## 1190 -74.75746 39.06786 -3
## 1191 -74.75746 39.06786 -3
## 1192 -74.75746 39.06786 -3
## 1193 -74.75746 39.06786 -3
## 1194 -74.75746 39.06786 -3
## 1195 -74.75746 39.06786 -3
## 1196 -74.75746 39.06786 -3
## 1197 -74.75746 39.06786 -3
## 1198 -74.75746 39.06786 -3
## 1199 -74.75746 39.06786 -3
## 1200 -74.75746 39.06786 -3
## 1201 -74.75746 39.06786 -3
## 1202 -74.75746 39.06786 -3
## 1203 -74.75746 39.06786 -3
## 1204 -74.75746 39.06786 -3
## 1205 -74.75746 39.06786 -3
## 1206 -74.75746 39.06786 -3
## 1207 -74.75746 39.06786 -3
## 1208 -74.75746 39.06786 -3
## 1209 -74.75746 39.06786 -3
## 1210 -74.75746 39.06786 -3
## 1211 -74.75746 39.06786 -3
## 1212 -74.75746 39.06786 -3
## 1213 -74.75746 39.06786 -3
## 1214 -74.75746 39.06786 -3
## 1215 -74.75746 39.06786 -3
## 1216 -74.75746 39.06786 -3
## 1217 -74.75746 39.06786 -3
## 1218 -74.75746 39.06786 -3
## 1219 -74.75746 39.06786 -3
## 1220 -74.75746 39.06786 -3
## 1221 -74.75746 39.06786 -3
## 1222 -74.75746 39.06786 -3
## 1223 -74.75746 39.06786 -3
## 1224 -74.75746 39.06786 -3
## 1225 -74.75746 39.06786 -3
## 1226 -74.75746 39.06786 -3
## 1227 -74.75746 39.06786 -3
## 1228 -74.75746 39.06786 -3
## 1229 -74.75746 39.06786 -3
## 1230 -74.75746 39.06786 -3
## 1231 -74.75746 39.06786 -3
## 1232 -74.75746 39.06786 -3
## 1233 -74.75746 39.06786 -3
## 1234 -74.75746 39.06786 -3
## 1235 -74.75746 39.06786 -3
## 1236 -70.97979 41.53583 -1
## 1237 -70.97979 41.53583 -1
## 1238 -70.97979 41.53583 -1
## 1239 -70.97979 41.53583 -1
## 1240 -70.97979 41.53583 -1
## 1241 -70.97979 41.53583 -1
## 1242 -70.97979 41.53583 -1
## 1243 -70.97979 41.53583 -1
## 1244 -70.97979 41.53583 -1
## 1245 -70.97979 41.53583 -1
## 1246 -70.97979 41.53583 -1
## 1247 -70.97979 41.53583 -1
## 1248 -70.97979 41.53583 -1
## 1249 -70.97979 41.53583 -1
## 1250 -70.97979 41.53583 -1
## 1251 -70.97979 41.53583 -1
## 1252 -70.97979 41.53583 -1
## 1253 -70.97979 41.53583 -1
## 1254 -70.97979 41.53583 -1
## 1255 -70.97979 41.53583 -1
## 1256 -70.97979 41.53583 -1
## 1257 -70.97979 41.53583 -1
## 1258 -70.97979 41.53583 -1
## 1259 -70.97979 41.53583 -1
## 1260 -70.97979 41.53583 -1
## 1261 -70.97979 41.53583 -1
## 1262 -70.97979 41.53583 -1
## 1263 -70.97979 41.53583 -1
## 1264 -70.97979 41.53583 -1
## 1265 -70.97979 41.53583 -1
## 1266 -71.52557 41.38235 -6
## 1267 -71.52557 41.38235 -6
## 1268 -71.52557 41.38235 -6
## 1269 -71.52557 41.38235 -6
## 1270 -71.52557 41.38235 -6
## 1271 -71.52557 41.38235 -6
## 1272 -71.52557 41.38235 -6
## 1273 -71.52557 41.38235 -6
## 1274 -71.52557 41.38235 -6
## 1275 -71.52557 41.38235 -6
## 1276 -71.52557 41.38235 -6
## 1277 -71.52557 41.38235 -6
## 1278 -71.52557 41.38235 -6
## 1279 -71.52557 41.38235 -6
## 1280 -71.52557 41.38235 -6
## 1281 -71.52557 41.38235 -6
## 1282 -71.52557 41.38235 -6
## 1283 -71.52557 41.38235 -6
## 1284 -71.52557 41.38235 -6
## 1285 -71.52557 41.38235 -6
## 1286 -71.52557 41.38235 -6
## 1287 -71.52557 41.38235 -6
## 1288 -71.52557 41.38235 -6
## 1289 -71.52557 41.38235 -6
## 1290 -71.52557 41.38235 -6
## 1291 -71.52557 41.38235 -6
## 1292 -71.52557 41.38235 -6
## 1293 -71.52557 41.38235 -6
## 1294 -71.52557 41.38235 -6
## 1295 -71.52557 41.38235 -6
## 1296 -71.52557 41.38235 -6
## 1297 -71.52557 41.38235 -6
## 1298 -71.52557 41.38235 -6
## 1299 -71.52557 41.38235 -6
## 1300 -71.52557 41.38235 -6
## 1301 -71.52557 41.38235 -6
## 1302 -71.52557 41.38235 -6
## 1303 -71.52557 41.38235 -6
## 1304 -71.52557 41.38235 -6
## 1305 -71.52557 41.38235 -6
## 1306 -71.52557 41.38235 -6
## 1307 -71.52557 41.38235 -6
## 1308 -71.52557 41.38235 -6
## 1309 -71.52557 41.38235 -6
## 1310 -71.52557 41.38235 -6
## 1311 -71.52557 41.38235 -6
## 1312 -71.52557 41.38235 -6
## 1313 -71.52557 41.38235 -6
## 1314 -70.90760 41.62486 -6
## 1315 -70.90760 41.62486 -6
## 1316 -70.90760 41.62486 -6
## 1317 -70.90760 41.62486 -6
## 1318 -70.90760 41.62486 -6
## 1319 -70.90760 41.62486 -6
## 1320 -70.90760 41.62486 -6
## 1321 -70.90760 41.62486 -6
## 1322 -70.90760 41.62486 -6
## 1323 -70.90760 41.62486 -6
## 1324 -70.90760 41.62486 -6
## 1325 -70.90760 41.62486 -6
## 1326 -70.90760 41.62486 -6
## 1327 -70.90760 41.62486 -6
## 1328 -70.90760 41.62486 -6
## 1329 -70.90760 41.62486 -6
## 1330 -70.90760 41.62486 -6
## 1331 -70.90760 41.62486 -6
## 1332 -70.90760 41.62486 -6
## 1333 -70.90760 41.62486 -6
## 1334 -70.90760 41.62486 -6
## 1335 -70.90760 41.62486 -6
## 1336 -70.90760 41.62486 -6
## 1337 -70.90760 41.62486 -6
## 1338 -70.90760 41.62486 -6
## 1339 -70.90760 41.62486 -6
## 1340 -70.90760 41.62486 -6
## 1341 -70.90760 41.62486 -6
## 1342 -70.90760 41.62486 -6
## 1343 -70.90760 41.62486 -6
## 1344 -77.92477 34.02835 0
## 1345 -77.92477 34.02835 0
## 1346 -77.92477 34.02835 0
## 1347 -77.92477 34.02835 0
## 1348 -77.92477 34.02835 0
## 1349 -77.92477 34.02835 0
## 1350 -77.92477 34.02835 0
## 1351 -77.92477 34.02835 0
## 1352 -77.92477 34.02835 0
## 1353 -77.92477 34.02835 0
## 1354 -77.92477 34.02835 0
## 1355 -77.92477 34.02835 0
## 1356 -77.92477 34.02835 0
## 1357 -77.92477 34.02835 0
## 1358 -77.92477 34.02835 0
## 1359 -77.92477 34.02835 0
## 1360 -77.92477 34.02835 0
## 1361 -77.92477 34.02835 0
## 1362 -77.92477 34.02835 0
## 1363 -77.92477 34.02835 0
## 1364 -77.92477 34.02835 0
## 1365 -77.92477 34.02835 0
## 1366 -77.92477 34.02835 0
## 1367 -77.92477 34.02835 0
## 1368 -77.92477 34.02835 0
## 1369 -77.92477 34.02835 0
## 1370 -77.92477 34.02835 0
## 1371 -77.92477 34.02835 0
## Warning in lc.dist(trans, ind_sites[, c(4, 3)], res = "dist"): One or more points are located outside of the [1;-15] depth range. You will get unrealistically huge distances unless you either increase the range of possible depths in trans.mat() or you move the problematic points in a spot where their depths fall within the [1;-15] depth range.
## You can use get.depth() to determine the depth of any point on a bathymetric map
#calculate dbmem and plot
dbmem_pos<-dbmem(distance_matrix, MEM.autocor = "positive")
## Warning in dudi.pco(matdist, scannf = FALSE, nf = 2): Non euclidean
## distance
plot(attributes(dbmem_pos)$values)
abline(h = 0)
The first 8 dbmems have vaules above 0 (i.e describe spatial autocorrelation)
temperatures
The environmental data is daily SSTs (NOAA OI SST V2 High Resolution Dataset) from
We extracted 10 years of daily optimum interpolated sea surface temperatures (OISSTs) from NOAA/OAR/ESRL PSD (Reynolds et al. 2007) for each of the sampling locations and calculated three derived environmental variables: mean temperature, mean annual minimum and mean annual maximum (table 5.1, fig. 5.1). The first principal component in a principal component analysis captured 90.1% of the variance among these variables, and all three variables were highly correlated as indicated by the biplot (supplementary fig. 5.1). While mean temperature and mean annual maximum demonstrate a linear relationship with latitude (R2 = 0.87 and 0.97; slope = -0.97 and - 128 1.1°C per ° latitude, respectively), minimum temperature is constrained by the freezing temperature of seawater and plateaued beyond 39 degrees of latitude. The value of the first temperature principal component and position along the shore for the sampling locations are also strongly correlated (R2 = 0.94, p < 10-9 ANOVA)
#env variable for pops
env<-read.table("./environmental/temp_summary.txt", sep = "\t", header = T)
env <- env %>%
left_join(., pops, by = "pop")
## Warning: Column `pop` joining factor and character vector, coercing into
## character vector
We use the covariance matrix calculated by PCANGSD to conduct a PCA that describes the genetic variation among samples. Then choose the number of PCs to retain using the Kaiser Guttman criteria
#pca
#to load pca: attach("./pcangsd/pca.r")
#snp_pca<-read.table("./pcangsd/pca.txt", header = T, sep= "\t")
#kaiser guttman
cutoff<-mean(pc1$sdev^2)
kg <- length((pc1$sdev^2)[(pc1$sdev^2)>cutoff])
#pcs to keep
snp_pcs <- pc1$x[,c(1:103)]
We keep the first kg pcs to fit our RDA on
Finally we get to running the RDA. We conduct a global RDA, an RDA for the temperature once conditioned on spatial data, and an RDA on the genetic DATA + spatial data once conditioned on the temperature. then we do model selection, significance testing and variance partitioning
#global rda
rda_null <- rda(snp_pcs ~ 1)
rda_ind<-rda(snp_pcs ~ as.matrix(env[,c(2)]) +as.matrix(env[,c(3)]) + as.matrix(env[,4]) + (as.matrix(dbmem_pos[,1])) + (as.matrix(dbmem_pos[,2])) + (as.matrix(dbmem_pos[,3])) + (as.matrix(dbmem_pos[,4]))+ (as.matrix(dbmem_pos[,5])) + (as.matrix(dbmem_pos[,6])) + (as.matrix(dbmem_pos[,7])) + (as.matrix(dbmem_pos[,8])))
#variable selection
ordi <- ordistep(rda_null, scope = formula(rda_ind))
#ordi$anova
#variance partitioning
vp <- varpart(snp_pcs ,as.matrix(env[,c(2,3,4)]) , (as.matrix(dbmem_pos)[,c(1:8)]))
#controlling for spatial autocor
rda_ind_con_spatial<-rda(snp_pcs ~ as.matrix(env[,c(2)]) + Condition((as.matrix(dbmem_pos)[,c(1:8)])))
#controlling for temp
rda_ind_con_temp<-rda(snp_pcs ~ Condition(as.matrix(env[,c(2)])) + (as.matrix(dbmem_pos)[,c(1:8)]))
#anovas
aov_global <- anova(rda_ind)
aov_con_spatial <- anova(rda_ind_con_spatial)
aov_con_temp <- anova(rda_ind_con_temp)
RDA Results Summary
This is a preliminary summary, more detailed results are to come.
We used the ordiR2step function to select our model. Using forward selection the full model was chosen (8 dbMEMs describing spatial autocorrelation and 3 environmental summary variables (mean temperature mean annual minimum, and mean annual maximum)).
Three RDAs were fit and tested for significance:
* global RDA - both the effect of spatial autocorrelation (the dbMEMs) and the effect of temerpature - p < 0.001
* RDA of temperature effect, conditioned on spatial autocorrelation - p < 0.001
* RDA of spatial effect conditioned on temperature - p < 0.001
Collectively these suggest that there is a significant relationship between the environmental temperature variable and the genetic covariance among samples, once the genetic data is conditioned on the effect of spatial autocorrelation AND vice versa. 76% of the total variation in the genetic dataset (or at least the first 103 of 1371 principal components of it) was constrained by the RDA.
Using the variance partitioning the relative effect of spatial autocorrelation and temperature could be estimated. Adjusted R2 for the effect of temperature, once spatial autocorrelation is accounted for, is 0.4%. In contrast 62% of the genetic dataset was explained by spatial autocorrelation alone. 14% of variation was explained jointly and could not be parsed. 23% was residual.
Examining the triplot suggests that the first dbmem separates northern from southern at the NJ.
RDA figures and results tables
#triplot
rda_sum<-summary(rda_ind)
scores<-as.data.frame(scores(rda_ind)$sites)
scores$pop<-pops$pop
arrows<-as.data.frame(rda_sum$biplot)
#order the factors
scores$pop <- factor( as.character(scores$pop), levels=c("GA", "SC", "WNC", "NC", "KC", "SH", "RB", "OC", "MG", "CT", "SR", "BP", "HB", "NBH", "P", "SLC", "ME", "MDIBL") )
ggplot()+geom_point(aes(RDA1, RDA2, color = pop), data = scores)+geom_segment(data=arrows, aes(x=0, y=0, xend=RDA1, yend=RDA2), arrow=arrow(length=unit(0.2,"cm")), alpha=0.75, color="red")+geom_text(data=arrows, aes(x=RDA1+0.2, y=RDA2, label=c("Max Temp", "Min Temp", "Mean Temp", "dbMEM1", "dbMEM2","dbMEM3","dbMEM4","dbMEM5", "dbMEM6","dbMEM7", "dbMEM8")), size = 5, vjust=1, color="red")+theme_bw()
rda_ind
## Call: rda(formula = snp_pcs ~ as.matrix(env[, c(2)]) +
## as.matrix(env[, c(3)]) + as.matrix(env[, 4]) +
## (as.matrix(dbmem_pos[, 1])) + (as.matrix(dbmem_pos[, 2])) +
## (as.matrix(dbmem_pos[, 3])) + (as.matrix(dbmem_pos[, 4])) +
## (as.matrix(dbmem_pos[, 5])) + (as.matrix(dbmem_pos[, 6])) +
## (as.matrix(dbmem_pos[, 7])) + (as.matrix(dbmem_pos[, 8])))
##
## Inertia Proportion Rank
## Total 6.4115 1.0000
## Constrained 4.9251 0.7682 11
## Unconstrained 1.4864 0.2318 103
## Inertia is variance
##
## Eigenvalues for constrained axes:
## RDA1 RDA2 RDA3 RDA4 RDA5 RDA6 RDA7 RDA8 RDA9 RDA10 RDA11
## 3.758 1.076 0.056 0.016 0.011 0.006 0.002 0.000 0.000 0.000 0.000
##
## Eigenvalues for unconstrained axes:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8
## 0.11497 0.05451 0.02401 0.02194 0.02174 0.02118 0.02082 0.02073
## (Showing 8 of 103 unconstrained eigenvalues)
aov_global
rda_ind_con_spatial
## Call: rda(formula = snp_pcs ~ as.matrix(env[, c(2)]) +
## Condition((as.matrix(dbmem_pos)[, c(1:8)])))
##
## Inertia Proportion Rank
## Total 6.4115155 1.0000000
## Conditional 4.8953346 0.7635222 8
## Constrained 0.0043957 0.0006856 1
## Unconstrained 1.5117852 0.2357922 103
## Inertia is variance
##
## Eigenvalues for constrained axes:
## RDA1
## 0.004396
##
## Eigenvalues for unconstrained axes:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8
## 0.13284 0.05644 0.02721 0.02201 0.02176 0.02119 0.02086 0.02074
## (Showing 8 of 103 unconstrained eigenvalues)
aov_con_spatial
rda_ind_con_temp
## Call: rda(formula = snp_pcs ~ Condition(as.matrix(env[, c(2)])) +
## (as.matrix(dbmem_pos)[, c(1:8)]))
##
## Inertia Proportion Rank
## Total 6.4115 1.0000
## Conditional 0.1808 0.0282 1
## Constrained 4.7189 0.7360 8
## Unconstrained 1.5118 0.2358 103
## Inertia is variance
##
## Eigenvalues for constrained axes:
## RDA1 RDA2 RDA3 RDA4 RDA5 RDA6 RDA7 RDA8
## 3.582 1.056 0.054 0.014 0.009 0.004 0.001 0.000
##
## Eigenvalues for unconstrained axes:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8
## 0.13284 0.05644 0.02721 0.02201 0.02176 0.02119 0.02086 0.02074
## (Showing 8 of 103 unconstrained eigenvalues)
aov_con_temp
notes for going forward
* it might be innappropirate to conduct anovas of rda condition on the combined environmental or dbmem matrices, it may be neccessary to indivudally test each variable, consitioned on all other variables… read into this before proceeding
* if we conduct the RDA on all pcs (not just those exceeding the kaiser guttman criteria) we can relate the RDA results more clearly to FST, i.e. we can see how much of the total variation among the cline is explained by RDA and compare this to results from amova and so on